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