source/wham/src/parmread.F

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