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