source/wham/src-M/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       call card_concat(controlcard,.true.)
  26       call reada(controlcard,'SCAL14',scal14,0.4d0)
  27       call reada(controlcard,'SCALSCP',scalscp,1.0d0)
  28       call reada(controlcard,'CUTOFF',cutoff_corr,7.0d0)
  29       call reada(controlcard,'DELT_CORR',delt_corr,0.5d0)
  30       r0_corr=cutoff_corr-delt_corr
  31       call readi(controlcard,"NRES",nres,0)
  32       iscode=index(controlcard,"ONE_LETTER")
  33       if (nres.le.0) then
  34         write (iout,*) "Error: no residues in molecule"
  35         return1
  36       endif
  37       if (nres.gt.maxres) then
  38         write (iout,*) "Error: too many residues",nres,maxres
  39       endif
  40       write(iout,*) 'nres=',nres
  41 C Read sequence of the protein
  42       if (iscode.gt.0) then
  43         read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
  44       else
  45         read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
  46       endif
  47 C Convert sequence to numeric code
  48       do i=1,nres
  49         itype(i)=rescode(i,sequence(i),iscode)
  50       enddo
  51       write (iout,*) "Numeric code:"
  52       write (iout,'(20i4)') (itype(i),i=1,nres)
  53       do i=1,nres-1
  54 #ifdef PROCOR
  55         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
  56 #else
  57         if (itype(i).eq.ntyp1) then
  58 #endif
  59           itel(i)=0
  60 #ifdef PROCOR
  61         else if (iabs(itype(i+1)).ne.20) then
  62 #else
  63         else if (iabs(itype(i)).ne.20) then
  64 #endif
  65           itel(i)=1
  66         else
  67           itel(i)=2
  68         endif
  69       enddo
  70        write (iout,*) "ITEL"
  71        do i=1,nres-1
  72          write (iout,*) i,itype(i),itel(i)
  73        enddo
  74       call read_bridge
  75
  76       if (with_dihed_constr) then
  77
  78       read (inp,*) ndih_constr
  79       if (ndih_constr.gt.0) then
  80 C        read (inp,*) ftors
  81 C        write (iout,*) 'FTORS',ftors
  82         read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i),
  83      &   i=1,ndih_constr)
  84         write (iout,*)
  85      &   'There are',ndih_constr,' constraints on phi angles.'
  86         do i=1,ndih_constr
  87           write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i),
  88      &    ftors(i)
  89         enddo
  90         do i=1,ndih_constr
  91           phi0(i)=deg2rad*phi0(i)
  92           drange(i)=deg2rad*drange(i)
  93         enddo
  94       endif
  95
  96       endif
  97       if (with_theta_constr) then
  98 C with_theta_constr is keyword allowing for occurance of theta constrains
  99       read (inp,*) ntheta_constr
 100 C ntheta_constr is the number of theta constrains
 101       if (ntheta_constr.gt.0) then
 102 C        read (inp,*) ftors
 103         read (inp,*) (itheta_constr(i),theta_constr0(i),
 104      &  theta_drange(i),for_thet_constr(i),
 105      &  i=1,ntheta_constr)
 106 C the above code reads from 1 to ntheta_constr
 107 C itheta_constr(i) residue i for which is theta_constr
 108 C theta_constr0 the global minimum value
 109 C theta_drange is range for which there is no energy penalty
 110 C for_thet_constr is the force constant for quartic energy penalty
 111 C E=k*x**4
 112 C        if(me.eq.king.or..not.out1file)then
 113          write (iout,*)
 114      &   'There are',ntheta_constr,' constraints on phi angles.'
 115          do i=1,ntheta_constr
 116           write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i),
 117      &    theta_drange(i),
 118      &    for_thet_constr(i)
 119          enddo
 120 C        endif
 121         do i=1,ntheta_constr
 122           theta_constr0(i)=deg2rad*theta_constr0(i)
 123           theta_drange(i)=deg2rad*theta_drange(i)
 124         enddo
 125 C        if(me.eq.king.or..not.out1file)
 126 C     &   write (iout,*) 'FTORS',ftors
 127 C        do i=1,ntheta_constr
 128 C          ii = itheta_constr(i)
 129 C          thetabound(1,ii) = phi0(i)-drange(i)
 130 C          thetabound(2,ii) = phi0(i)+drange(i)
 131 C        enddo
 132       endif ! ntheta_constr.gt.0
 133       endif! with_theta_constr
 134       nnt=1
 135       nct=nres
 136       if (itype(1).eq.ntyp1) nnt=2
 137       if (itype(nres).eq.ntyp1) nct=nct-1
 138       write(iout,*) 'NNT=',NNT,' NCT=',NCT
 139       call setup_var
 140       call init_int_table
 141       if (ns.gt.0) then
 142         write (iout,'(/a,i3,a)') 'The chain contains',ns,
 143      &  ' disulfide-bridging cysteines.'
 144         write (iout,'(20i4)') (iss(i),i=1,ns)
 145        if (dyn_ss) then
 146           write(iout,*)"Running with dynamic disulfide-bond formation"
 147        else
 148         write (iout,'(/a/)') 'Pre-formed links are:'
 149         do i=1,nss
 150           i1=ihpb(i)-nres
 151           i2=jhpb(i)-nres
 152           it1=itype(i1)
 153           it2=itype(i2)
 154          write (iout,'(2a,i3,3a,i3,a,3f10.3)')
 155      &    restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',
 156      &    dhpb(i),ebr,forcon(i)
 157         enddo
 158       endif
 159       endif
 160       write (iout,'(a)')
 161       if (ns.gt.0.and.dyn_ss) then
 162           do i=nss+1,nhpb
 163             ihpb(i-nss)=ihpb(i)
 164             jhpb(i-nss)=jhpb(i)
 165             forcon(i-nss)=forcon(i)
 166             dhpb(i-nss)=dhpb(i)
 167           enddo
 168           nhpb=nhpb-nss
 169           nss=0
 170           call hpb_partition
 171           do i=1,ns
 172             dyn_ss_mask(iss(i))=.true.
 173           enddo
 174       endif
 175       return
 176       end
 177 c-----------------------------------------------------------------------------
 178       logical function seq_comp(itypea,itypeb,length)
 179       implicit none
 180       integer length,itypea(length),itypeb(length)
 181       integer i
 182       do i=1,length
 183         if (itypea(i).ne.itypeb(i)) then
 184           seq_comp=.false.
 185           return
 186         endif
 187       enddo
 188       seq_comp=.true.
 189       return
 190       end
 191 c-----------------------------------------------------------------------------
 192       subroutine read_bridge
 193 C Read information about disulfide bridges.
 194       implicit real*8 (a-h,o-z)
 195       include 'DIMENSIONS'
 196       include 'DIMENSIONS.ZSCOPT'
 197       include 'COMMON.IOUNITS'
 198       include 'COMMON.GEO'
 199       include 'COMMON.VAR'
 200       include 'COMMON.INTERACT'
 201       include 'COMMON.NAMES'
 202       include 'COMMON.CHAIN'
 203       include 'COMMON.FFIELD'
 204       include 'COMMON.SBRIDGE'
 205 C Read bridging residues.
 206       read (inp,*) ns,(iss(i),i=1,ns)
 207       print *,'ns=',ns
 208       write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
 209 C Check whether the specified bridging residues are cystines.
 210       do i=1,ns
 211         if (itype(iss(i)).ne.1) then
 212           write (iout,'(2a,i3,a)')
 213      &   'Do you REALLY think that the residue ',
 214      &    restyp(itype(iss(i))),i,
 215      &   ' can form a disulfide bridge?!!!'
 216           write (*,'(2a,i3,a)')
 217      &   'Do you REALLY think that the residue ',
 218      &    restyp(itype(iss(i))),i,
 219      &   ' can form a disulfide bridge?!!!'
 220          stop
 221         endif
 222       enddo
 223 C Read preformed bridges.
 224       if (ns.gt.0) then
 225       read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
 226       write (iout,*) 'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
 227       if (nss.gt.0) then
 228         nhpb=nss
 229 C Check if the residues involved in bridges are in the specified list of
 230 C bridging residues.
 231         do i=1,nss
 232           do j=1,i-1
 233             if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
 234      &      .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
 235               write (iout,'(a,i3,a)') 'Disulfide pair',i,
 236      &      ' contains residues present in other pairs.'
 237               write (*,'(a,i3,a)') 'Disulfide pair',i,
 238      &      ' contains residues present in other pairs.'
 239               stop
 240             endif
 241           enddo
 242           do j=1,ns
 243             if (ihpb(i).eq.iss(j)) goto 10
 244           enddo
 245           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
 246    10     continue
 247           do j=1,ns
 248             if (jhpb(i).eq.iss(j)) goto 20
 249           enddo
 250           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
 251    20     continue
 252           dhpb(i)=dbr
 253           forcon(i)=fbr
 254         enddo
 255         do i=1,nss
 256           ihpb(i)=ihpb(i)+nres
 257           jhpb(i)=jhpb(i)+nres
 258         enddo
 259       endif
 260       endif
 261       return
 262       end
 263 c------------------------------------------------------------------------------
 264       subroutine read_angles(kanal,iscor,energ,iprot,*)
 265       implicit real*8 (a-h,o-z)
 266       include 'DIMENSIONS'
 267       include 'DIMENSIONS.ZSCOPT'
 268       include 'COMMON.INTERACT'
 269       include 'COMMON.SBRIDGE'
 270       include 'COMMON.GEO'
 271       include 'COMMON.VAR'
 272       include 'COMMON.CHAIN'
 273       include 'COMMON.IOUNITS'
 274       character*80 lineh
 275       read(kanal,'(a80)',end=10,err=10) lineh
 276       read(lineh(:5),*,err=8) ic
 277       read(lineh(6:),*,err=8) energ
 278       goto 9
 279     8 ic=1
 280       print *,'error, assuming e=1d10',lineh
 281       energ=1d10
 282       nss=0
 283     9 continue
 284       read(lineh(18:),*,end=10,err=10) nss
 285       IF (NSS.LT.9) THEN
 286         read (lineh(20:),*,end=10,err=10)
 287      &  (IHPB(I),JHPB(I),I=1,NSS),iscor
 288       ELSE
 289         read (lineh(20:),*,end=10,err=10) (IHPB(I),JHPB(I),I=1,8)
 290         read (kanal,*,end=10,err=10) (IHPB(I),JHPB(I),
 291      &    I=9,NSS),iscor
 292       ENDIF
 293 c      print *,"energy",energ," iscor",iscor
 294       read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
 295       read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
 296       read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
 297       read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
 298       do i=1,nres
 299         theta(i)=deg2rad*theta(i)
 300         phi(i)=deg2rad*phi(i)
 301         alph(i)=deg2rad*alph(i)
 302         omeg(i)=deg2rad*omeg(i)
 303       enddo
 304       return
 305    10 return1
 306       end