source/cluster/wham/src-M/parmread.F

   1       subroutine parmread
   2 C
   3 C Read the parameters of the probability distributions of the virtual-bond
   4 C valence angles and the side chains and energy parameters.
   5 C
   6       implicit real*8 (a-h,o-z)
   7       include 'DIMENSIONS'
   8       include 'sizesclu.dat'
   9       include 'COMMON.IOUNITS'
  10       include 'COMMON.CHAIN'
  11       include 'COMMON.INTERACT'
  12       include 'COMMON.GEO'
  13       include 'COMMON.LOCAL'
  14       include 'COMMON.TORSION'
  15       include 'COMMON.FFIELD'
  16       include 'COMMON.NAMES'
  17       include 'COMMON.SBRIDGE'
  18       include 'COMMON.SCCOR'
  19       include 'COMMON.SCROT'
  20       character*1 t1,t2,t3
  21       character*1 onelett(-2:2) /"p","a","G","A","P"/
  22       logical lprint
  23       dimension blower(3,3,maxlob)
  24       double precision ip,mp
  25 C
  26 C Body
  27 C
  28 C Set LPRINT=.TRUE. for debugging
  29       dwa16=2.0d0**(1.0d0/6.0d0)
  30       lprint=.false.
  31       itypro=20
  32 C Assign virtual-bond length
  33       vbl=3.8D0
  34       vblinv=1.0D0/vbl
  35       vblinv2=vblinv*vblinv
  36 #ifdef CRYST_BOND
  37       read (ibond,*) vbldp0,akp
  38       do i=1,ntyp
  39         nbondterm(i)=1
  40         read (ibond,*) vbldsc0(1,i),aksc(1,i)
  41         dsc(i) = vbldsc0(1,i)
  42         if (i.eq.10) then
  43           dsc_inv(i)=0.0D0
  44         else
  45           dsc_inv(i)=1.0D0/dsc(i)
  46         endif
  47       enddo
  48 #else
  49       read (ibond,*) ijunk,vbldp0,akp,rjunk
  50       do i=1,ntyp
  51         read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),
  52      &   j=1,nbondterm(i))
  53         dsc(i) = vbldsc0(1,i)
  54         if (i.eq.10) then
  55           dsc_inv(i)=0.0D0
  56         else
  57           dsc_inv(i)=1.0D0/dsc(i)
  58         endif
  59       enddo
  60 #endif
  61       if (lprint) then
  62         write(iout,'(/a/)')"Force constants virtual bonds:"
  63         write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K',
  64      &   'inertia','Pstok'
  65         write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0
  66         do i=1,ntyp
  67           write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),
  68      &      vbldsc0(1,i),aksc(1,i),abond0(1,i)
  69           do j=2,nbondterm(i)
  70             write (iout,'(13x,3f10.5)')
  71      &        vbldsc0(j,i),aksc(j,i),abond0(j,i)
  72           enddo
  73         enddo
  74       endif
  75 #ifdef CRYST_THETA
  76 C
  77 C Read the parameters of the probability distribution/energy expression
  78 C of the virtual-bond valence angles theta
  79 C
  80       do i=1,ntyp
  81         read (ithep,*) a0thet(i),(athet(j,i,1,1),j=1,2),
  82      &    (bthet(j,i,1,1),j=1,2)
  83         read (ithep,*) (polthet(j,i),j=0,3)
  84         read (ithep,*) (gthet(j,i),j=1,3)
  85         read (ithep,*) theta0(i),sig0(i),sigc0(i)
  86         sigc0(i)=sigc0(i)**2
  87       enddo
  88       do i=1,ntyp
  89       athet(1,i,1,-1)=athet(1,i,1,1)
  90       athet(2,i,1,-1)=athet(2,i,1,1)
  91       bthet(1,i,1,-1)=-bthet(1,i,1,1)
  92       bthet(2,i,1,-1)=-bthet(2,i,1,1)
  93       athet(1,i,-1,1)=-athet(1,i,1,1)
  94       athet(2,i,-1,1)=-athet(2,i,1,1)
  95       bthet(1,i,-1,1)=bthet(1,i,1,1)
  96       bthet(2,i,-1,1)=bthet(2,i,1,1)
  97       enddo
  98       do i=-ntyp,-1
  99       a0thet(i)=a0thet(-i)
 100       athet(1,i,-1,-1)=athet(1,-i,1,1)
 101       athet(2,i,-1,-1)=-athet(2,-i,1,1)
 102       bthet(1,i,-1,-1)=bthet(1,-i,1,1)
 103       bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
 104       athet(1,i,-1,1)=athet(1,-i,1,1)
 105       athet(2,i,-1,1)=-athet(2,-i,1,1)
 106       bthet(1,i,-1,1)=-bthet(1,-i,1,1)
 107       bthet(2,i,-1,1)=bthet(2,-i,1,1)
 108       athet(1,i,1,-1)=-athet(1,-i,1,1)
 109       athet(2,i,1,-1)=athet(2,-i,1,1)
 110       bthet(1,i,1,-1)=bthet(1,-i,1,1)
 111       bthet(2,i,1,-1)=-bthet(2,-i,1,1)
 112       theta0(i)=theta0(-i)
 113       sig0(i)=sig0(-i)
 114       sigc0(i)=sigc0(-i)
 115        do j=0,3
 116         polthet(j,i)=polthet(j,-i)
 117        enddo
 118        do j=1,3
 119          gthet(j,i)=gthet(j,-i)
 120        enddo
 121       enddo
 122       close (ithep)
 123       if (lprint) then
 124 c       write (iout,'(a)')
 125 c    &    'Parameters of the virtual-bond valence angles:'
 126 c       write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
 127 c    & '    ATHETA0   ','         A1   ','        A2    ',
 128 c    & '        B1    ','         B2   '
 129 c       do i=1,ntyp
 130 c         write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
 131 c    &        a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2)
 132 c       enddo
 133 c       write (iout,'(/a/9x,5a/79(1h-))')
 134 c    & 'Parameters of the expression for sigma(theta_c):',
 135 c    & '     ALPH0    ','      ALPH1   ','     ALPH2    ',
 136 c    & '     ALPH3    ','    SIGMA0C   '
 137 c       do i=1,ntyp
 138 c         write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
 139 c    &      (polthet(j,i),j=0,3),sigc0(i)
 140 c       enddo
 141 c       write (iout,'(/a/9x,5a/79(1h-))')
 142 c    & 'Parameters of the second gaussian:',
 143 c    & '    THETA0    ','     SIGMA0   ','        G1    ',
 144 c    & '        G2    ','         G3   '
 145 c       do i=1,ntyp
 146 c         write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
 147 c    &       sig0(i),(gthet(j,i),j=1,3)
 148 c       enddo
 149         write (iout,'(a)')
 150      &    'Parameters of the virtual-bond valence angles:'
 151         write (iout,'(/a/9x,5a/79(1h-))')
 152      & 'Coefficients of expansion',
 153      & '     theta0   ','    a1*10^2   ','   a2*10^2    ',
 154      & '   b1*10^1    ','    b2*10^1   '
 155         do i=1,ntyp
 156           write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),
 157      &        a0thet(i),(100*athet(j,i,1,1),j=1,2),
 158      &        (10*bthet(j,i,1,1),j=1,2)
 159         enddo
 160         write (iout,'(/a/9x,5a/79(1h-))')
 161      & 'Parameters of the expression for sigma(theta_c):',
 162      & ' alpha0       ','  alph1       ',' alph2        ',
 163      & ' alhp3        ','   sigma0c    '
 164         do i=1,ntyp
 165           write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),
 166      &      (polthet(j,i),j=0,3),sigc0(i)
 167         enddo
 168         write (iout,'(/a/9x,5a/79(1h-))')
 169      & 'Parameters of the second gaussian:',
 170      & '    theta0    ','  sigma0*10^2 ','      G1*10^-1',
 171      & '        G2    ','   G3*10^1    '
 172         do i=1,ntyp
 173           write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),
 174      &       100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
 175         enddo
 176       endif
 177 #else
 178 C
 179 C Read the parameters of Utheta determined from ab initio surfaces
 180 C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
 181 C
 182       read (ithep,*) nthetyp,ntheterm,ntheterm2,
 183      &  ntheterm3,nsingle,ndouble
 184       nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
 185       read (ithep,*) (ithetyp(i),i=1,ntyp1)
 186       do i=1,maxthetyp
 187         do j=1,maxthetyp
 188           do k=1,maxthetyp
 189             aa0thet(i,j,k)=0.0d0
 190             do l=1,ntheterm
 191               aathet(l,i,j,k)=0.0d0
 192             enddo
 193             do l=1,ntheterm2
 194               do m=1,nsingle
 195                 bbthet(m,l,i,j,k)=0.0d0
 196                 ccthet(m,l,i,j,k)=0.0d0
 197                 ddthet(m,l,i,j,k)=0.0d0
 198                 eethet(m,l,i,j,k)=0.0d0
 199               enddo
 200             enddo
 201             do l=1,ntheterm3
 202               do m=1,ndouble
 203                 do mm=1,ndouble
 204                  ffthet(mm,m,l,i,j,k)=0.0d0
 205                  ggthet(mm,m,l,i,j,k)=0.0d0
 206                 enddo
 207               enddo
 208             enddo
 209           enddo
 210         enddo
 211       enddo
 212       do i=1,nthetyp
 213         do j=1,nthetyp
 214           do k=1,nthetyp
 215             read (ithep,'(3a)') res1,res2,res3
 216             read (ithep,*) aa0thet(i,j,k)
 217             read (ithep,*)(aathet(l,i,j,k),l=1,ntheterm)
 218             read (ithep,*)
 219      &       ((bbthet(lll,ll,i,j,k),lll=1,nsingle),
 220      &        (ccthet(lll,ll,i,j,k),lll=1,nsingle),
 221      &        (ddthet(lll,ll,i,j,k),lll=1,nsingle),
 222      &        (eethet(lll,ll,i,j,k),lll=1,nsingle),ll=1,ntheterm2)
 223             read (ithep,*)
 224      &      (((ffthet(llll,lll,ll,i,j,k),ffthet(lll,llll,ll,i,j,k),
 225      &         ggthet(llll,lll,ll,i,j,k),ggthet(lll,llll,ll,i,j,k),
 226      &         llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
 227           enddo
 228         enddo
 229       enddo
 230 C
 231 C For dummy ends assign glycine-type coefficients of theta-only terms; the
 232 C coefficients of theta-and-gamma-dependent terms are zero.
 233 C
 234       do i=1,nthetyp
 235         do j=1,nthetyp
 236           do l=1,ntheterm
 237             aathet(l,i,j,nthetyp+1)=aathet(l,i,j,1)
 238             aathet(l,nthetyp+1,i,j)=aathet(l,1,i,j)
 239           enddo
 240           aa0thet(i,j,nthetyp+1)=aa0thet(i,j,1)
 241           aa0thet(nthetyp+1,i,j)=aa0thet(1,i,j)
 242         enddo
 243         do l=1,ntheterm
 244           aathet(l,nthetyp+1,i,nthetyp+1)=aathet(l,1,i,1)
 245         enddo
 246         aa0thet(nthetyp+1,i,nthetyp+1)=aa0thet(1,i,1)
 247       enddo
 248 C
 249 C Control printout of the coefficients of virtual-bond-angle potentials
 250 C
 251       if (lprint) then
 252         write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
 253         do i=1,nthetyp+1
 254           do j=1,nthetyp+1
 255             do k=1,nthetyp+1
 256               write (iout,'(//4a)')
 257      &         'Type ',onelett(i),onelett(j),onelett(k)
 258               write (iout,'(//a,10x,a)') " l","a[l]"
 259               write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k)
 260               write (iout,'(i2,1pe15.5)')
 261      &           (l,aathet(l,i,j,k),l=1,ntheterm)
 262             do l=1,ntheterm2
 263               write (iout,'(//2h m,4(9x,a,3h[m,i1,1h]))')
 264      &          "b",l,"c",l,"d",l,"e",l
 265               do m=1,nsingle
 266                 write (iout,'(i2,4(1pe15.5))') m,
 267      &          bbthet(m,l,i,j,k),ccthet(m,l,i,j,k),
 268      &          ddthet(m,l,i,j,k),eethet(m,l,i,j,k)
 269               enddo
 270             enddo
 271             do l=1,ntheterm3
 272               write (iout,'(//3hm,n,4(6x,a,5h[m,n,i1,1h]))')
 273      &          "f+",l,"f-",l,"g+",l,"g-",l
 274               do m=2,ndouble
 275                 do n=1,m-1
 276                   write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
 277      &              ffthet(n,m,l,i,j,k),ffthet(m,n,l,i,j,k),
 278      &              ggthet(n,m,l,i,j,k),ggthet(m,n,l,i,j,k)
 279                 enddo
 280               enddo
 281             enddo
 282           enddo
 283         enddo
 284       enddo
 285       call flush(iout)
 286       endif
 287 #endif
 288
 289 #ifdef CRYST_SC
 290 C
 291 C Read the parameters of the probability distribution/energy expression
 292 C of the side chains.
 293 C
 294       do i=1,ntyp
 295         read (irotam,'(3x,i3,f8.3)') nlob(i),dsc(i)
 296         if (i.eq.10) then
 297           dsc_inv(i)=0.0D0
 298         else
 299           dsc_inv(i)=1.0D0/dsc(i)
 300         endif
 301         if (i.ne.10) then
 302         do j=1,nlob(i)
 303           do k=1,3
 304             do l=1,3
 305               blower(l,k,j)=0.0D0
 306             enddo
 307           enddo
 308         enddo
 309         bsc(1,i)=0.0D0
 310         read(irotam,*)(censc(k,1,i),k=1,3),((blower(k,l,1),l=1,k),k=1,3)
 311         censc(1,1,-i)=censc(1,1,i)
 312         censc(2,1,-i)=censc(2,1,i)
 313         censc(3,1,-i)=-censc(3,1,i)
 314         do j=2,nlob(i)
 315           read (irotam,*) bsc(j,i)
 316           read (irotam,*) (censc(k,j,i),k=1,3),
 317      &                                 ((blower(k,l,j),l=1,k),k=1,3)
 318         censc(1,j,-i)=censc(1,j,i)
 319         censc(2,j,-i)=censc(2,j,i)
 320         censc(3,j,-i)=-censc(3,j,i)
 321 C BSC is amplitude of Gaussian
 322         enddo
 323         do j=1,nlob(i)
 324           do k=1,3
 325             do l=1,k
 326               akl=0.0D0
 327               do m=1,3
 328                 akl=akl+blower(k,m,j)*blower(l,m,j)
 329               enddo
 330               gaussc(k,l,j,i)=akl
 331               gaussc(l,k,j,i)=akl
 332              if (((k.eq.3).and.(l.ne.3))
 333      &        .or.((l.eq.3).and.(k.ne.3))) then
 334                 gaussc(k,l,j,-i)=-akl
 335                 gaussc(l,k,j,-i)=-akl
 336               else
 337                 gaussc(k,l,j,-i)=akl
 338                 gaussc(l,k,j,-i)=akl
 339               endif
 340             enddo
 341           enddo
 342         enddo
 343         endif
 344       enddo
 345       close (irotam)
 346       if (lprint) then
 347         write (iout,'(/a)') 'Parameters of side-chain local geometry'
 348         do i=1,ntyp
 349           nlobi=nlob(i)
 350           if (nlobi.gt.0) then
 351           write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
 352      &     ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
 353 c          write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
 354 c          write (iout,'(a,f10.4,4(16x,f10.4))')
 355 c     &                             'Center  ',(bsc(j,i),j=1,nlobi)
 356 c          write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),j=1,nlobi)
 357            write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
 358      &                             'log h',(bsc(j,i),j=1,nlobi)
 359            write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
 360      &    'x',((censc(k,j,i),k=1,3),j=1,nlobi)
 361 c          write (iout,'(a)')
 362 c         do j=1,nlobi
 363 c           ind=0
 364 c           do k=1,3
 365 c             do l=1,k
 366 c              ind=ind+1
 367 c              blower(k,l,j)=gaussc(ind,j,i)
 368 c             enddo
 369 c           enddo
 370 c         enddo
 371           do k=1,3
 372             write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
 373      &                 ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
 374           enddo
 375           endif
 376         enddo
 377       endif
 378 #else
 379 C
 380 C Read scrot parameters for potentials determined from all-atom AM1 calculations
 381 C added by Urszula Kozlowska 07/11/2007
 382 C
 383       do i=1,ntyp
 384         read (irotam,*)
 385        if (i.eq.10) then
 386          read (irotam,*)
 387        else
 388          do j=1,65
 389            read(irotam,*) sc_parmin(j,i)
 390          enddo
 391        endif
 392       enddo
 393 #endif
 394       close(irotam)
 395 #ifdef CRYST_TOR
 396 C
 397 C Read torsional parameters in old format
 398 C
 399       read (itorp,*) ntortyp,nterm_old
 400       write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
 401       read (itorp,*) (itortyp(i),i=1,ntyp)
 402       do i=1,ntortyp
 403         do j=1,ntortyp
 404           read (itorp,'(a)')
 405           do k=1,nterm_old
 406             read (itorp,*) kk,v1(k,j,i),v2(k,j,i)
 407           enddo
 408         enddo
 409       enddo
 410       close (itorp)
 411       if (lprint) then
 412         write (iout,'(/a/)') 'Torsional constants:'
 413         do i=1,ntortyp
 414           do j=1,ntortyp
 415             write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
 416             write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
 417           enddo
 418         enddo
 419       endif
 420 #else
 421 C
 422 C Read torsional parameters
 423 C
 424       read (itorp,*) ntortyp
 425       read (itorp,*) (itortyp(i),i=1,ntyp)
 426       write (iout,*) 'ntortyp',ntortyp
 427       do iblock=1,2
 428       do i=-ntyp,-1
 429        itortyp(i)=-itortyp(-i)
 430       enddo
 431 c      write (iout,*) 'ntortyp',ntortyp
 432       do i=0,ntortyp-1
 433         do j=-ntortyp+1,ntortyp-1
 434           read (itorp,*) nterm(i,j,iblock),
 435      &          nlor(i,j,iblock)
 436           nterm(-i,-j,iblock)=nterm(i,j,iblock)
 437           nlor(-i,-j,iblock)=nlor(i,j,iblock)
 438           v0ij=0.0d0
 439           si=-1.0d0
 440           do k=1,nterm(i,j,iblock)
 441             read (itorp,*) kk,v1(k,i,j,iblock),
 442      &      v2(k,i,j,iblock)
 443             v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
 444             v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
 445             v0ij=v0ij+si*v1(k,i,j,iblock)
 446             si=-si
 447           enddo
 448           do k=1,nlor(i,j,iblock)
 449             read (itorp,*) kk,vlor1(k,i,j),
 450      &        vlor2(k,i,j),vlor3(k,i,j)
 451             v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
 452           enddo
 453           v0(i,j,iblock)=v0ij
 454           v0(-i,-j,iblock)=v0ij
 455         enddo
 456       enddo
 457       enddo
 458       close (itorp)
 459       if (lprint) then
 460         write (iout,'(/a/)') 'Torsional constants:'
 461         do i=1,ntortyp
 462           do j=1,ntortyp
 463             write (iout,*) 'ityp',i,' jtyp',j
 464             write (iout,*) 'Fourier constants'
 465             do k=1,nterm(i,j,iblock)
 466              write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
 467      &        v2(k,i,j,iblock)
 468             enddo
 469             write (iout,*) 'Lorenz constants'
 470             do k=1,nlor(i,j,iblock)
 471               write (iout,'(3(1pe15.5))')
 472      &         vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
 473             enddo
 474           enddo
 475         enddo
 476       endif
 477 C
 478 C 6/23/01 Read parameters for double torsionals
 479 C
 480       do iblock=1,2
 481       do i=0,ntortyp-1
 482         do j=-ntortyp+1,ntortyp-1
 483           do k=-ntortyp+1,ntortyp-1
 484             read (itordp,'(3a1)') t1,t2,t3
 485             if (t1.ne.onelett(i) .or. t2.ne.onelett(j)
 486      &        .or. t3.ne.onelett(k)) then
 487               write (iout,*) "Error in double torsional parameter file",
 488      &         i,j,k,t1,t2,t3
 489                stop "Error in double torsional parameter file"
 490             endif
 491           read (itordp,*) ntermd_1(i,j,k,iblock),
 492      &         ntermd_2(i,j,k,iblock)
 493             ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
 494             ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
 495             read (itordp,*) (v1c(1,l,i,j,k,iblock),l=1,
 496      &         ntermd_1(i,j,k,iblock))
 497             read (itordp,*) (v1s(1,l,i,j,k,iblock),l=1,
 498      &         ntermd_1(i,j,k,iblock))
 499             read (itordp,*) (v1c(2,l,i,j,k,iblock),l=1,
 500      &         ntermd_1(i,j,k,iblock))
 501             read (itordp,*) (v1s(2,l,i,j,k,iblock),l=1,
 502      &         ntermd_1(i,j,k,iblock))
 503 C Martix of D parameters for one dimesional foureir series
 504             do l=1,ntermd_1(i,j,k,iblock)
 505              v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
 506              v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
 507              v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
 508              v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
 509 c            write(iout,*) "whcodze" ,
 510 c     & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
 511             enddo
 512             read (itordp,*) ((v2c(l,m,i,j,k,iblock),
 513      &         v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
 514      &         v2s(m,l,i,j,k,iblock),
 515      &         m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
 516 C Martix of D parameters for two dimesional fourier series
 517             do l=1,ntermd_2(i,j,k,iblock)
 518              do m=1,l-1
 519              v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
 520              v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
 521              v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
 522              v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
 523              enddo!m
 524             enddo!l
 525           enddo!k
 526         enddo!j
 527       enddo!i
 528       enddo!iblock
 529       if (lprint) then
 530       write (iout,*)
 531       write (iout,*) 'Constants for double torsionals'
 532       do iblock=1,2
 533       do i=0,ntortyp-1
 534         do j=-ntortyp+1,ntortyp-1
 535           do k=-ntortyp+1,ntortyp-1
 536             write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
 537      &        ' nsingle',ntermd_1(i,j,k,iblock),
 538      &        ' ndouble',ntermd_2(i,j,k,iblock)
 539             write (iout,*)
 540             write (iout,*) 'Single angles:'
 541             do l=1,ntermd_1(i,j,k,iblock)
 542               write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
 543      &           v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
 544      &           v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),
 545      &           v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
 546             enddo
 547             write (iout,*)
 548             write (iout,*) 'Pairs of angles:'
 549             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
 550             do l=1,ntermd_2(i,j,k,iblock)
 551               write (iout,'(i5,20f10.5)')
 552      &         l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
 553             enddo
 554             write (iout,*)
 555            write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
 556             do l=1,ntermd_2(i,j,k,iblock)
 557               write (iout,'(i5,20f10.5)')
 558      &         l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),
 559      &         (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
 560             enddo
 561             write (iout,*)
 562           enddo
 563         enddo
 564       enddo
 565       enddo
 566       endif
 567 #endif
 568 C
 569 C 5/21/07 (AL) Read coefficients of the backbone-local sidechain-local
 570 C         correlation energies.
 571 C
 572       read (isccor,*) nterm_sccor
 573       do i=1,20
 574         do j=1,20
 575           read (isccor,'(a)')
 576           do k=1,nterm_sccor
 577             read (isccor,*)
 578      &        kk,v1sccor(k,i,j),v2sccor(k,i,j)
 579           enddo
 580         enddo
 581       enddo
 582       close (isccor)
 583       if (lprint) then
 584         write (iout,'(/a/)') 'Torsional constants of SCCORR:'
 585         do i=1,20
 586           do j=1,20
 587             write (iout,*) 'ityp',i,' jtyp',j
 588             do k=1,nterm_sccor
 589               write (iout,'(2(1pe15.5))') v1sccor(k,i,j),v2sccor(k,i,j)
 590             enddo
 591           enddo
 592         enddo
 593       endif
 594 C
 595 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
 596 C         interaction energy of the Gly, Ala, and Pro prototypes.
 597 C
 598       read (ifourier,*) nloctyp
 599       do i=0,nloctyp-1
 600         read (ifourier,*)
 601         read (ifourier,*) (b(ii,i),ii=1,13)
 602         if (lprint) then
 603         write (iout,*) 'Type',i
 604         write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
 605         endif
 606         B1(1,i)  = b(3,i)
 607         B1(2,i)  = b(5,i)
 608         B1(1,-i) = b(3,i)
 609         B1(2,-i) = -b(5,i)
 610         B1tilde(1,i) = b(3,i)
 611         B1tilde(2,i) =-b(5,i)
 612         B1tilde(1,-i) =-b(3,i)
 613         B1tilde(2,-i) =b(5,i)
 614         B2(1,i)  = b(2,i)
 615         B2(2,i)  = b(4,i)
 616         B2(1,-i)  =b(2,i)
 617         B2(2,-i)  =-b(4,i)
 618         CC(1,1,i)= b(7,i)
 619         CC(2,2,i)=-b(7,i)
 620         CC(2,1,i)= b(9,i)
 621         CC(1,2,i)= b(9,i)
 622         CC(1,1,-i)= b(7,i)
 623         CC(2,2,-i)=-b(7,i)
 624         CC(2,1,-i)=-b(9,i)
 625         CC(1,2,-i)=-b(9,i)
 626         Ctilde(1,1,i)=b(7,i)
 627         Ctilde(1,2,i)=b(9,i)
 628         Ctilde(2,1,i)=-b(9,i)
 629         Ctilde(2,2,i)=b(7,i)
 630         Ctilde(1,1,-i)=b(7,i)
 631         Ctilde(1,2,-i)=-b(9,i)
 632         Ctilde(2,1,-i)=b(9,i)
 633         Ctilde(2,2,-i)=b(7,i)
 634         DD(1,1,i)= b(6,i)
 635         DD(2,2,i)=-b(6,i)
 636         DD(2,1,i)= b(8,i)
 637         DD(1,2,i)= b(8,i)
 638         DD(1,1,-i)= b(6,i)
 639         DD(2,2,-i)=-b(6,i)
 640         DD(2,1,-i)=-b(8,i)
 641         DD(1,2,-i)=-b(8,i)
 642         Dtilde(1,1,i)=b(6,i)
 643         Dtilde(1,2,i)=b(8,i)
 644         Dtilde(2,1,i)=-b(8,i)
 645         Dtilde(2,2,i)=b(6,i)
 646         Dtilde(1,1,-i)=b(6,i)
 647         Dtilde(1,2,-i)=-b(8,i)
 648         Dtilde(2,1,-i)=b(8,i)
 649         Dtilde(2,2,-i)=b(6,i)
 650         EE(1,1,i)= b(10,i)+b(11,i)
 651         EE(2,2,i)=-b(10,i)+b(11,i)
 652         EE(2,1,i)= b(12,i)-b(13,i)
 653         EE(1,2,i)= b(12,i)+b(13,i)
 654         EE(1,1,-i)= b(10,i)+b(11,i)
 655         EE(2,2,-i)=-b(10,i)+b(11,i)
 656         EE(2,1,-i)=-b(12,i)+b(13,i)
 657         EE(1,2,-i)=-b(12,i)-b(13,i)
 658       enddo
 659       if (lprint) then
 660       do i=1,nloctyp
 661         write (iout,*) 'Type',i
 662         write (iout,*) 'B1'
 663 c        write (iout,'(f10.5)') B1(:,i)
 664         write(iout,*) B1(1,i),B1(2,i)
 665         write (iout,*) 'B2'
 666 c        write (iout,'(f10.5)') B2(:,i)
 667         write(iout,*) B2(1,i),B2(2,i)
 668         write (iout,*) 'CC'
 669         do j=1,2
 670           write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
 671         enddo
 672         write(iout,*) 'DD'
 673         do j=1,2
 674           write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
 675         enddo
 676         write(iout,*) 'EE'
 677         do j=1,2
 678           write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
 679         enddo
 680       enddo
 681       endif
 682 C
 683 C Read electrostatic-interaction parameters
 684 C
 685       if (lprint) then
 686         write (iout,'(/a)') 'Electrostatic interaction constants:'
 687         write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
 688      &            'IT','JT','APP','BPP','AEL6','AEL3'
 689       endif
 690       read (ielep,*) ((epp(i,j),j=1,2),i=1,2)
 691       read (ielep,*) ((rpp(i,j),j=1,2),i=1,2)
 692       read (ielep,*) ((elpp6(i,j),j=1,2),i=1,2)
 693       read (ielep,*) ((elpp3(i,j),j=1,2),i=1,2)
 694       close (ielep)
 695       do i=1,2
 696         do j=1,2
 697         rri=rpp(i,j)**6
 698         app (i,j)=epp(i,j)*rri*rri
 699         bpp (i,j)=-2.0D0*epp(i,j)*rri
 700         ael6(i,j)=elpp6(i,j)*4.2D0**6
 701         ael3(i,j)=elpp3(i,j)*4.2D0**3
 702         if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
 703      &                    ael6(i,j),ael3(i,j)
 704         enddo
 705       enddo
 706 C
 707 C Read side-chain interaction parameters.
 708 C
 709       read (isidep,*) ipot,expon
 710       if (ipot.lt.1 .or. ipot.gt.5) then
 711         write (iout,'(2a)') 'Error while reading SC interaction',
 712      &               'potential file - unknown potential type.'
 713         stop
 714       endif
 715       expon2=expon/2
 716       write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
 717      & ', exponents are ',expon,2*expon
 718       goto (10,20,30,30,40) ipot
 719 C----------------------- LJ potential ---------------------------------
 720    10 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
 721       if (lprint) then
 722         write (iout,'(/a/)') 'Parameters of the LJ potential:'
 723         write (iout,'(a/)') 'The epsilon array:'
 724         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 725         write (iout,'(/a)') 'One-body parameters:'
 726         write (iout,'(a,4x,a)') 'residue','sigma'
 727         write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
 728       endif
 729       goto 50
 730 C----------------------- LJK potential --------------------------------
 731    20 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
 732      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
 733       if (lprint) then
 734         write (iout,'(/a/)') 'Parameters of the LJK potential:'
 735         write (iout,'(a/)') 'The epsilon array:'
 736         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 737         write (iout,'(/a)') 'One-body parameters:'
 738         write (iout,'(a,4x,2a)') 'residue','   sigma  ','    r0    '
 739         write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
 740      &        i=1,ntyp)
 741       endif
 742       goto 50
 743 C---------------------- GB or BP potential -----------------------------
 744    30 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
 745      &  (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip0(i),i=1,ntyp),
 746      &  (alp(i),i=1,ntyp)
 747 C For the GB potential convert sigma'**2 into chi'
 748       if (ipot.eq.4) then
 749         do i=1,ntyp
 750           chip(i)=(chip0(i)-1.0D0)/(chip0(i)+1.0D0)
 751         enddo
 752       endif
 753       if (lprint) then
 754         write (iout,'(/a/)') 'Parameters of the BP potential:'
 755         write (iout,'(a/)') 'The epsilon array:'
 756         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 757         write (iout,'(/a)') 'One-body parameters:'
 758         write (iout,'(a,4x,4a)') 'residue','   sigma  ','s||/s_|_^2',
 759      &       '    chip  ','    alph  '
 760         write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
 761      &                     chip(i),alp(i),i=1,ntyp)
 762       endif
 763       goto 50
 764 C--------------------- GBV potential -----------------------------------
 765    40 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
 766      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
 767      &  (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
 768       if (lprint) then
 769         write (iout,'(/a/)') 'Parameters of the GBV potential:'
 770         write (iout,'(a/)') 'The epsilon array:'
 771         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 772         write (iout,'(/a)') 'One-body parameters:'
 773         write (iout,'(a,4x,5a)') 'residue','   sigma  ','    r0    ',
 774      &      's||/s_|_^2','    chip  ','    alph  '
 775         write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
 776      &           sigii(i),chip(i),alp(i),i=1,ntyp)
 777       endif
 778    50 continue
 779       close (isidep)
 780 C-----------------------------------------------------------------------
 781 C Calculate the "working" parameters of SC interactions.
 782       do i=2,ntyp
 783         do j=1,i-1
 784           eps(i,j)=eps(j,i)
 785         enddo
 786       enddo
 787       do i=1,ntyp
 788         do j=i,ntyp
 789           sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
 790           sigma(j,i)=sigma(i,j)
 791           rs0(i,j)=dwa16*sigma(i,j)
 792           rs0(j,i)=rs0(i,j)
 793         enddo
 794       enddo
 795       if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
 796      & 'Working parameters of the SC interactions:',
 797      & '     a    ','     b    ','   augm   ','  sigma ','   r0   ',
 798      & '  chi1   ','   chi2   '
 799       do i=1,ntyp
 800         do j=i,ntyp
 801           epsij=eps(i,j)
 802           if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
 803             rrij=sigma(i,j)
 804           else
 805             rrij=rr0(i)+rr0(j)
 806           endif
 807           r0(i,j)=rrij
 808           r0(j,i)=rrij
 809           rrij=rrij**expon
 810           epsij=eps(i,j)
 811           sigeps=dsign(1.0D0,epsij)
 812           epsij=dabs(epsij)
 813           aa(i,j)=epsij*rrij*rrij
 814           bb(i,j)=-sigeps*epsij*rrij
 815           aa(j,i)=aa(i,j)
 816           bb(j,i)=bb(i,j)
 817           if (ipot.gt.2) then
 818             sigt1sq=sigma0(i)**2
 819             sigt2sq=sigma0(j)**2
 820             sigii1=sigii(i)
 821             sigii2=sigii(j)
 822             ratsig1=sigt2sq/sigt1sq
 823             ratsig2=1.0D0/ratsig1
 824             chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
 825             if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
 826             rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
 827           else
 828             rsum_max=sigma(i,j)
 829           endif
 830 c         if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
 831             sigmaii(i,j)=rsum_max
 832             sigmaii(j,i)=rsum_max
 833 c         else
 834 c           sigmaii(i,j)=r0(i,j)
 835 c           sigmaii(j,i)=r0(i,j)
 836 c         endif
 837 cd        write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
 838           if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
 839             r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
 840             augm(i,j)=epsij*r_augm**(2*expon)
 841 c           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
 842             augm(j,i)=augm(i,j)
 843           else
 844             augm(i,j)=0.0D0
 845             augm(j,i)=0.0D0
 846           endif
 847           if (lprint) then
 848             write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
 849      &      restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
 850      &      sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
 851           endif
 852         enddo
 853       enddo
 854 C
 855 C Define the SC-p interaction constants
 856 C
 857       do i=1,20
 858         do j=1,2
 859           eps_scp(i,j)=-1.5d0
 860           rscp(i,j)=4.0d0
 861         enddo
 862       enddo
 863 #ifdef OLDSCP
 864       do i=1,20
 865 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
 866 C helix formation)
 867 c       aad(i,1)=0.3D0*4.0D0**12
 868 C Following line for constants currently implemented
 869 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
 870         aad(i,1)=1.5D0*4.0D0**12
 871 c       aad(i,1)=0.17D0*5.6D0**12
 872         aad(i,2)=aad(i,1)
 873 C "Soft" SC-p repulsion
 874         bad(i,1)=0.0D0
 875 C Following line for constants currently implemented
 876 c       aad(i,1)=0.3D0*4.0D0**6
 877 C "Hard" SC-p repulsion
 878         bad(i,1)=3.0D0*4.0D0**6
 879 c       bad(i,1)=-2.0D0*0.17D0*5.6D0**6
 880         bad(i,2)=bad(i,1)
 881 c       aad(i,1)=0.0D0
 882 c       aad(i,2)=0.0D0
 883 c       bad(i,1)=1228.8D0
 884 c       bad(i,2)=1228.8D0
 885       enddo
 886 #else
 887 C
 888 C 8/9/01 Read the SC-p interaction constants from file
 889 C
 890       do i=1,ntyp
 891         read (iscpp,*) (eps_scp(i,j),rscp(i,j),j=1,2)
 892       enddo
 893       do i=1,ntyp
 894         aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
 895         aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
 896         bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
 897         bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
 898       enddo
 899
 900       if (lprint) then
 901         write (iout,*) "Parameters of SC-p interactions:"
 902         do i=1,20
 903           write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
 904      &     eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
 905         enddo
 906       endif
 907 #endif
 908 C
 909 C Define the constants of the disulfide bridge
 910 C
 911       ebr=-5.50D0
 912 c
 913 c Old arbitrary potential - commented out.
 914 c
 915 c      dbr= 4.20D0
 916 c      fbr= 3.30D0
 917 c
 918 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
 919 c energy surface of diethyl disulfide.
 920 c A. Liwo and U. Kozlowska, 11/24/03
 921 c
 922       D0CM = 3.78d0
 923       AKCM = 15.1d0
 924       AKTH = 11.0d0
 925       AKCT = 12.0d0
 926       V1SS =-1.08d0
 927       V2SS = 7.61d0
 928       V3SS = 13.7d0
 929
 930       write (iout,'(/a)') "Disulfide bridge parameters:"
 931       write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
 932       write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
 933       write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
 934       write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
 935      & ' v3ss:',v3ss
 936       return
 937       end