source/cluster/wham/src/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 Read of Side-chain backbone correlation parameters
 569 C Modified 11 May 2012 by Adasko
 570 CCC
 571 C
 572       read (isccor,*) nsccortyp
 573       read (isccor,*) (isccortyp(i),i=1,ntyp)
 574 c      write (iout,*) 'ntortyp',ntortyp
 575       do i=-ntyp,-1
 576         isccortyp(i)=-isccortyp(-i)
 577       enddo
 578       iscprol=isccortyp(20)
 579       maxinter=3
 580 cc maxinter is maximum interaction sites
 581       do l=1,maxinter
 582       do i=1,nsccortyp
 583         do j=1,nsccortyp
 584           read (isccor,*)
 585      &nterm_sccor(i,j),nlor_sccor(i,j)
 586           v0ijsccor=0.0d0
 587           v0ijsccor1=0.0d0
 588           v0ijsccor2=0.0d0
 589           v0ijsccor3=0.0d0
 590          si=-1.0d0
 591           nterm_sccor(-i,j)=nterm_sccor(i,j)
 592           nterm_sccor(-i,-j)=nterm_sccor(i,j)
 593           nterm_sccor(i,-j)=nterm_sccor(i,j)
 594           do k=1,nterm_sccor(i,j)
 595             read (isccor,*) kk,v1sccor(k,l,i,j)
 596      &    ,v2sccor(k,l,i,j)
 597             if (j.eq.iscprol) then
 598              if (i.eq.isccortyp(10)) then
 599              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
 600              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
 601              else
 602              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
 603      &                        +v2sccor(k,l,i,j)*dsqrt(0.75d0)
 604             v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
 605      &                        +v1sccor(k,l,i,j)*dsqrt(0.75d0)
 606               v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
 607               v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
 608               v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
 609               v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
 610               endif
 611              else
 612               if (i.eq.isccortyp(10)) then
 613               v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
 614               v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
 615               else
 616                if (j.eq.isccortyp(10)) then
 617              v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
 618              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
 619                else
 620              v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
 621              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
 622              v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
 623              v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
 624              v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
 625              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
 626                 endif
 627                endif
 628             endif
 629             v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
 630             v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
 631             v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
 632             v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
 633             si=-si
 634           enddo
 635           do k=1,nlor_sccor(i,j)
 636             read (isccor,*) kk,vlor1sccor(k,i,j),
 637      &        vlor2sccor(k,i,j),vlor3sccor(k,i,j)
 638             v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
 639      &(1+vlor3sccor(k,i,j)**2)
 640           enddo
 641            v0sccor(l,i,j)=v0ijsccor
 642            v0sccor(l,-i,j)=v0ijsccor1
 643            v0sccor(l,i,-j)=v0ijsccor2
 644            v0sccor(l,-i,-j)=v0ijsccor3
 645         enddo
 646       enddo
 647       enddo
 648       close (isccor)
 649
 650       if (lprint) then
 651         write (iout,'(/a/)') 'Torsional constants:'
 652         do i=1,nsccortyp
 653           do j=1,nsccortyp
 654             write (iout,*) 'ityp',i,' jtyp',j
 655             write (iout,*) 'Fourier constants'
 656             do k=1,nterm_sccor(i,j)
 657               write (iout,'(2(1pe15.5))')
 658      & v1sccor(k,l,i,j),v2sccor(k,l,i,j)
 659             enddo
 660             write (iout,*) 'Lorenz constants'
 661             do k=1,nlor_sccor(i,j)
 662               write (iout,'(3(1pe15.5))')
 663      &         vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
 664             enddo
 665           enddo
 666         enddo
 667       endif
 668
 669 C
 670 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
 671 C         interaction energy of the Gly, Ala, and Pro prototypes.
 672 C
 673       read (ifourier,*) nloctyp
 674       do i=0,nloctyp-1
 675         read (ifourier,*)
 676         read (ifourier,*) (b(ii,i),ii=1,13)
 677         if (lprint) then
 678         write (iout,*) 'Type',i
 679         write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
 680         endif
 681         B1(1,i)  = b(3,i)
 682         B1(2,i)  = b(5,i)
 683         B1(1,-i) = b(3,i)
 684         B1(2,-i) = -b(5,i)
 685         B1tilde(1,i) = b(3,i)
 686         B1tilde(2,i) =-b(5,i)
 687         B1tilde(1,-i) =-b(3,i)
 688         B1tilde(2,-i) =b(5,i)
 689         B2(1,i)  = b(2,i)
 690         B2(2,i)  = b(4,i)
 691         B2(1,-i)  =b(2,i)
 692         B2(2,-i)  =-b(4,i)
 693         CC(1,1,i)= b(7,i)
 694         CC(2,2,i)=-b(7,i)
 695         CC(2,1,i)= b(9,i)
 696         CC(1,2,i)= b(9,i)
 697         CC(1,1,-i)= b(7,i)
 698         CC(2,2,-i)=-b(7,i)
 699         CC(2,1,-i)=-b(9,i)
 700         CC(1,2,-i)=-b(9,i)
 701         Ctilde(1,1,i)=b(7,i)
 702         Ctilde(1,2,i)=b(9,i)
 703         Ctilde(2,1,i)=-b(9,i)
 704         Ctilde(2,2,i)=b(7,i)
 705         Ctilde(1,1,-i)=b(7,i)
 706         Ctilde(1,2,-i)=-b(9,i)
 707         Ctilde(2,1,-i)=b(9,i)
 708         Ctilde(2,2,-i)=b(7,i)
 709         DD(1,1,i)= b(6,i)
 710         DD(2,2,i)=-b(6,i)
 711         DD(2,1,i)= b(8,i)
 712         DD(1,2,i)= b(8,i)
 713         DD(1,1,-i)= b(6,i)
 714         DD(2,2,-i)=-b(6,i)
 715         DD(2,1,-i)=-b(8,i)
 716         DD(1,2,-i)=-b(8,i)
 717         Dtilde(1,1,i)=b(6,i)
 718         Dtilde(1,2,i)=b(8,i)
 719         Dtilde(2,1,i)=-b(8,i)
 720         Dtilde(2,2,i)=b(6,i)
 721         Dtilde(1,1,-i)=b(6,i)
 722         Dtilde(1,2,-i)=-b(8,i)
 723         Dtilde(2,1,-i)=b(8,i)
 724         Dtilde(2,2,-i)=b(6,i)
 725         EE(1,1,i)= b(10,i)+b(11,i)
 726         EE(2,2,i)=-b(10,i)+b(11,i)
 727         EE(2,1,i)= b(12,i)-b(13,i)
 728         EE(1,2,i)= b(12,i)+b(13,i)
 729         EE(1,1,-i)= b(10,i)+b(11,i)
 730         EE(2,2,-i)=-b(10,i)+b(11,i)
 731         EE(2,1,-i)=-b(12,i)+b(13,i)
 732         EE(1,2,-i)=-b(12,i)-b(13,i)
 733       enddo
 734       if (lprint) then
 735       do i=1,nloctyp
 736         write (iout,*) 'Type',i
 737         write (iout,*) 'B1'
 738 c        write (iout,'(f10.5)') B1(:,i)
 739         write(iout,*) B1(1,i),B1(2,i)
 740         write (iout,*) 'B2'
 741 c        write (iout,'(f10.5)') B2(:,i)
 742         write(iout,*) B2(1,i),B2(2,i)
 743         write (iout,*) 'CC'
 744         do j=1,2
 745           write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
 746         enddo
 747         write(iout,*) 'DD'
 748         do j=1,2
 749           write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
 750         enddo
 751         write(iout,*) 'EE'
 752         do j=1,2
 753           write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
 754         enddo
 755       enddo
 756       endif
 757 C
 758 C Read electrostatic-interaction parameters
 759 C
 760       if (lprint) then
 761         write (iout,'(/a)') 'Electrostatic interaction constants:'
 762         write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
 763      &            'IT','JT','APP','BPP','AEL6','AEL3'
 764       endif
 765       read (ielep,*) ((epp(i,j),j=1,2),i=1,2)
 766       read (ielep,*) ((rpp(i,j),j=1,2),i=1,2)
 767       read (ielep,*) ((elpp6(i,j),j=1,2),i=1,2)
 768       read (ielep,*) ((elpp3(i,j),j=1,2),i=1,2)
 769       close (ielep)
 770       do i=1,2
 771         do j=1,2
 772         rri=rpp(i,j)**6
 773         app (i,j)=epp(i,j)*rri*rri
 774         bpp (i,j)=-2.0D0*epp(i,j)*rri
 775         ael6(i,j)=elpp6(i,j)*4.2D0**6
 776         ael3(i,j)=elpp3(i,j)*4.2D0**3
 777         if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
 778      &                    ael6(i,j),ael3(i,j)
 779         enddo
 780       enddo
 781 C
 782 C Read side-chain interaction parameters.
 783 C
 784       read (isidep,*) ipot,expon
 785       if (ipot.lt.1 .or. ipot.gt.5) then
 786         write (iout,'(2a)') 'Error while reading SC interaction',
 787      &               'potential file - unknown potential type.'
 788         stop
 789       endif
 790       expon2=expon/2
 791       write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
 792      & ', exponents are ',expon,2*expon
 793       goto (10,20,30,30,40) ipot
 794 C----------------------- LJ potential ---------------------------------
 795    10 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
 796       if (lprint) then
 797         write (iout,'(/a/)') 'Parameters of the LJ potential:'
 798         write (iout,'(a/)') 'The epsilon array:'
 799         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 800         write (iout,'(/a)') 'One-body parameters:'
 801         write (iout,'(a,4x,a)') 'residue','sigma'
 802         write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
 803       endif
 804       goto 50
 805 C----------------------- LJK potential --------------------------------
 806    20 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
 807      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
 808       if (lprint) then
 809         write (iout,'(/a/)') 'Parameters of the LJK potential:'
 810         write (iout,'(a/)') 'The epsilon array:'
 811         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 812         write (iout,'(/a)') 'One-body parameters:'
 813         write (iout,'(a,4x,2a)') 'residue','   sigma  ','    r0    '
 814         write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
 815      &        i=1,ntyp)
 816       endif
 817       goto 50
 818 C---------------------- GB or BP potential -----------------------------
 819    30 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
 820      &  (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip0(i),i=1,ntyp),
 821      &  (alp(i),i=1,ntyp)
 822 C For the GB potential convert sigma'**2 into chi'
 823       if (ipot.eq.4) then
 824         do i=1,ntyp
 825           chip(i)=(chip0(i)-1.0D0)/(chip0(i)+1.0D0)
 826         enddo
 827       endif
 828       if (lprint) then
 829         write (iout,'(/a/)') 'Parameters of the BP potential:'
 830         write (iout,'(a/)') 'The epsilon array:'
 831         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 832         write (iout,'(/a)') 'One-body parameters:'
 833         write (iout,'(a,4x,4a)') 'residue','   sigma  ','s||/s_|_^2',
 834      &       '    chip  ','    alph  '
 835         write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
 836      &                     chip(i),alp(i),i=1,ntyp)
 837       endif
 838       goto 50
 839 C--------------------- GBV potential -----------------------------------
 840    40 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
 841      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
 842      &  (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
 843       if (lprint) then
 844         write (iout,'(/a/)') 'Parameters of the GBV potential:'
 845         write (iout,'(a/)') 'The epsilon array:'
 846         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 847         write (iout,'(/a)') 'One-body parameters:'
 848         write (iout,'(a,4x,5a)') 'residue','   sigma  ','    r0    ',
 849      &      's||/s_|_^2','    chip  ','    alph  '
 850         write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
 851      &           sigii(i),chip(i),alp(i),i=1,ntyp)
 852       endif
 853    50 continue
 854       close (isidep)
 855 C-----------------------------------------------------------------------
 856 C Calculate the "working" parameters of SC interactions.
 857       do i=2,ntyp
 858         do j=1,i-1
 859           eps(i,j)=eps(j,i)
 860         enddo
 861       enddo
 862       do i=1,ntyp
 863         do j=i,ntyp
 864           sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
 865           sigma(j,i)=sigma(i,j)
 866           rs0(i,j)=dwa16*sigma(i,j)
 867           rs0(j,i)=rs0(i,j)
 868         enddo
 869       enddo
 870       if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
 871      & 'Working parameters of the SC interactions:',
 872      & '     a    ','     b    ','   augm   ','  sigma ','   r0   ',
 873      & '  chi1   ','   chi2   '
 874       do i=1,ntyp
 875         do j=i,ntyp
 876           epsij=eps(i,j)
 877           if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
 878             rrij=sigma(i,j)
 879           else
 880             rrij=rr0(i)+rr0(j)
 881           endif
 882           r0(i,j)=rrij
 883           r0(j,i)=rrij
 884           rrij=rrij**expon
 885           epsij=eps(i,j)
 886           sigeps=dsign(1.0D0,epsij)
 887           epsij=dabs(epsij)
 888           aa(i,j)=epsij*rrij*rrij
 889           bb(i,j)=-sigeps*epsij*rrij
 890           aa(j,i)=aa(i,j)
 891           bb(j,i)=bb(i,j)
 892           if (ipot.gt.2) then
 893             sigt1sq=sigma0(i)**2
 894             sigt2sq=sigma0(j)**2
 895             sigii1=sigii(i)
 896             sigii2=sigii(j)
 897             ratsig1=sigt2sq/sigt1sq
 898             ratsig2=1.0D0/ratsig1
 899             chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
 900             if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
 901             rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
 902           else
 903             rsum_max=sigma(i,j)
 904           endif
 905 c         if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
 906             sigmaii(i,j)=rsum_max
 907             sigmaii(j,i)=rsum_max
 908 c         else
 909 c           sigmaii(i,j)=r0(i,j)
 910 c           sigmaii(j,i)=r0(i,j)
 911 c         endif
 912 cd        write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
 913           if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
 914             r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
 915             augm(i,j)=epsij*r_augm**(2*expon)
 916 c           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
 917             augm(j,i)=augm(i,j)
 918           else
 919             augm(i,j)=0.0D0
 920             augm(j,i)=0.0D0
 921           endif
 922           if (lprint) then
 923             write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
 924      &      restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
 925      &      sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
 926           endif
 927         enddo
 928       enddo
 929 C
 930 C Define the SC-p interaction constants
 931 C
 932       do i=1,20
 933         do j=1,2
 934           eps_scp(i,j)=-1.5d0
 935           rscp(i,j)=4.0d0
 936         enddo
 937       enddo
 938 #ifdef OLDSCP
 939       do i=1,20
 940 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
 941 C helix formation)
 942 c       aad(i,1)=0.3D0*4.0D0**12
 943 C Following line for constants currently implemented
 944 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
 945         aad(i,1)=1.5D0*4.0D0**12
 946 c       aad(i,1)=0.17D0*5.6D0**12
 947         aad(i,2)=aad(i,1)
 948 C "Soft" SC-p repulsion
 949         bad(i,1)=0.0D0
 950 C Following line for constants currently implemented
 951 c       aad(i,1)=0.3D0*4.0D0**6
 952 C "Hard" SC-p repulsion
 953         bad(i,1)=3.0D0*4.0D0**6
 954 c       bad(i,1)=-2.0D0*0.17D0*5.6D0**6
 955         bad(i,2)=bad(i,1)
 956 c       aad(i,1)=0.0D0
 957 c       aad(i,2)=0.0D0
 958 c       bad(i,1)=1228.8D0
 959 c       bad(i,2)=1228.8D0
 960       enddo
 961 #else
 962 C
 963 C 8/9/01 Read the SC-p interaction constants from file
 964 C
 965       do i=1,ntyp
 966         read (iscpp,*) (eps_scp(i,j),rscp(i,j),j=1,2)
 967       enddo
 968       do i=1,ntyp
 969         aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
 970         aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
 971         bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
 972         bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
 973       enddo
 974
 975       if (lprint) then
 976         write (iout,*) "Parameters of SC-p interactions:"
 977         do i=1,20
 978           write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
 979      &     eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
 980         enddo
 981       endif
 982 #endif
 983 C
 984 C Define the constants of the disulfide bridge
 985 C
 986       ebr=-5.50D0
 987 c
 988 c Old arbitrary potential - commented out.
 989 c
 990 c      dbr= 4.20D0
 991 c      fbr= 3.30D0
 992 c
 993 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
 994 c energy surface of diethyl disulfide.
 995 c A. Liwo and U. Kozlowska, 11/24/03
 996 c
 997       D0CM = 3.78d0
 998       AKCM = 15.1d0
 999       AKTH = 11.0d0
1000       AKCT = 12.0d0
1001       V1SS =-1.08d0
1002       V2SS = 7.61d0
1003       V3SS = 13.7d0
1004
1005       write (iout,'(/a)') "Disulfide bridge parameters:"
1006       write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1007       write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1008       write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1009       write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1010      & ' v3ss:',v3ss
1011       return
1012       end