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