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=-ntyp1,-1
 179         ithetyp(i)=-ithetyp(-i)
 180       enddo
 181       do iblock=1,2
 182       do i=-maxthetyp,maxthetyp
 183         do j=-maxthetyp,maxthetyp
 184           do k=-maxthetyp,maxthetyp
 185             aa0thet(i,j,k,iblock)=0.0d0
 186             do l=1,ntheterm
 187               aathet(l,i,j,k,iblock)=0.0d0
 188             enddo
 189             do l=1,ntheterm2
 190               do m=1,nsingle
 191                 bbthet(m,l,i,j,k,iblock)=0.0d0
 192                 ccthet(m,l,i,j,k,iblock)=0.0d0
 193                 ddthet(m,l,i,j,k,iblock)=0.0d0
 194                 eethet(m,l,i,j,k,iblock)=0.0d0
 195               enddo
 196             enddo
 197             do l=1,ntheterm3
 198               do m=1,ndouble
 199                 do mm=1,ndouble
 200                  ffthet(mm,m,l,i,j,k,iblock)=0.0d0
 201                  ggthet(mm,m,l,i,j,k,iblock)=0.0d0
 202                 enddo
 203               enddo
 204             enddo
 205           enddo
 206         enddo
 207       enddo
 208       enddo
 209 c VAR:iblock means terminally blocking group 1=non-proline 2=proline
 210       do iblock=1,2
 211 c VAR:ntethtyp is type of theta potentials type currently 0=glycine
 212 c VAR:1=non-glicyne non-proline 2=proline
 213 c VAR:negative values for D-aminoacid
 214       do i=0,nthetyp
 215         do j=-nthetyp,nthetyp
 216           do k=-nthetyp,nthetyp
 217             read (ithep,'(6a)',end=111,err=111) res1
 218 c VAR: aa0thet is variable describing the average value of Foureir
 219 c VAR: expansion series
 220             read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock)
 221 c VAR: aathet is foureir expansion in theta/2 angle for full formula
 222 c VAR: look at the fitting equation in Kozlowska et al., J. Phys.:
 223 Condens. Matter 19 (2007) 285203 and Sieradzan et al., unpublished
 224             read (ithep,*,end=111,err=111)
 225      &(aathet(l,i,j,k,iblock),l=1,ntheterm)
 226             read (ithep,*,end=111,err=111)
 227      &       ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
 228      &        (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
 229      &        (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
 230      &        (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle),
 231      &        ll=1,ntheterm2)
 232             read (ithep,*,end=111,err=111)
 233      &      (((ffthet(llll,lll,ll,i,j,k,iblock),
 234      &         ffthet(lll,llll,ll,i,j,k,iblock),
 235      &         ggthet(llll,lll,ll,i,j,k,iblock),
 236      &         ggthet(lll,llll,ll,i,j,k,iblock),
 237      &         llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
 238           enddo
 239         enddo
 240       enddo
 241
 242
 243 C
 244 C For dummy ends assign glycine-type coefficients of theta-only terms; the
 245 C coefficients of theta-and-gamma-dependent terms are zero.
 246 C IF YOU WANT VALENCE POTENTIALS FOR DUMMY ATOM UNCOMENT BELOW (NOT
 247 C RECOMENTDED AFTER VERSION 3.3)
 248 c      do i=1,nthetyp
 249 c        do j=1,nthetyp
 250 c          do l=1,ntheterm
 251 c            aathet(l,i,j,nthetyp+1,iblock)=aathet(l,i,j,1,iblock)
 252 c            aathet(l,nthetyp+1,i,j,iblock)=aathet(l,1,i,j,iblock)
 253 c          enddo
 254 c          aa0thet(i,j,nthetyp+1,iblock)=aa0thet(i,j,1,iblock)
 255 c          aa0thet(nthetyp+1,i,j,iblock)=aa0thet(1,i,j,iblock)
 256 c        enddo
 257 c        do l=1,ntheterm
 258 c          aathet(l,nthetyp+1,i,nthetyp+1,iblock)=aathet(l,1,i,1,iblock)
 259 c        enddo
 260 c        aa0thet(nthetyp+1,i,nthetyp+1,iblock)=aa0thet(1,i,1,iblock)
 261 c      enddo
 262 c      enddo
 263 C AND COMMENT THE LOOPS BELOW
 264       do i=1,nthetyp
 265         do j=1,nthetyp
 266           do l=1,ntheterm
 267             aathet(l,i,j,nthetyp+1,iblock)=0.0d0
 268             aathet(l,nthetyp+1,i,j,iblock)=0.0d0
 269           enddo
 270           aa0thet(i,j,nthetyp+1,iblock)=0.0d0
 271           aa0thet(nthetyp+1,i,j,iblock)=0.0d0
 272         enddo
 273         do l=1,ntheterm
 274           aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
 275         enddo
 276         aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
 277       enddo
 278       enddo
 279 C TILL HERE
 280 C Substitution for D aminoacids from symmetry.
 281       do iblock=1,2
 282       do i=-nthetyp,0
 283         do j=-nthetyp,nthetyp
 284           do k=-nthetyp,nthetyp
 285            aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
 286            do l=1,ntheterm
 287            aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
 288            enddo
 289            do ll=1,ntheterm2
 290             do lll=1,nsingle
 291             bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
 292             ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
 293             ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
 294             eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
 295             enddo
 296           enddo
 297           do ll=1,ntheterm3
 298            do lll=2,ndouble
 299             do llll=1,lll-1
 300             ffthet(llll,lll,ll,i,j,k,iblock)=
 301      &      ffthet(llll,lll,ll,-i,-j,-k,iblock)
 302             ffthet(lll,llll,ll,i,j,k,iblock)=
 303      &      ffthet(lll,llll,ll,-i,-j,-k,iblock)
 304             ggthet(llll,lll,ll,i,j,k,iblock)=
 305      &      -ggthet(llll,lll,ll,-i,-j,-k,iblock)
 306             ggthet(lll,llll,ll,i,j,k,iblock)=
 307      &      -ggthet(lll,llll,ll,-i,-j,-k,iblock)
 308             enddo !ll
 309            enddo  !lll
 310           enddo   !llll
 311          enddo    !k
 312         enddo     !j
 313        enddo      !i
 314       enddo       !iblock
 315 C
 316 C Control printout of the coefficients of virtual-bond-angle potentials
 317 C
 318       if (lprint) then
 319         write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
 320         do i=1,nthetyp+1
 321           do j=1,nthetyp+1
 322             do k=1,nthetyp+1
 323               write (iout,'(//4a)')
 324      &         'Type ',onelett(i),onelett(j),onelett(k)
 325               write (iout,'(//a,10x,a)') " l","a[l]"
 326               write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
 327               write (iout,'(i2,1pe15.5)')
 328      &           (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
 329             do l=1,ntheterm2
 330               write (iout,'(//2h m,4(9x,a,3h[m,,i1,1h]))')
 331      &          "b",l,"c",l,"d",l,"e",l
 332               do m=1,nsingle
 333                 write (iout,'(i2,4(1pe15.5))') m,
 334      &          bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),
 335      &          ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
 336               enddo
 337             enddo
 338             do l=1,ntheterm3
 339               write (iout,'(//3hm,n,4(6x,a,5h[m,n,,i1,1h]))')
 340      &          "f+",l,"f-",l,"g+",l,"g-",l
 341               do m=2,ndouble
 342                 do n=1,m-1
 343                   write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
 344      &              ffthet(n,m,l,i,j,k,iblock),
 345      &              ffthet(m,n,l,i,j,k,iblock),
 346      &              ggthet(n,m,l,i,j,k,iblock),
 347      &              ggthet(m,n,l,i,j,k,iblock)
 348                 enddo
 349               enddo
 350             enddo
 351           enddo
 352         enddo
 353       enddo
 354       call flush(iout)
 355       endif
 356       write (2,*) "Start reading THETA_PDB"
 357       do i=1,ntyp
 358         read (ithep_pdb,*,err=111,end=111)
 359      &     a0thet(i),(athet(j,i,1,1),j=1,2),
 360      &    (bthet(j,i,1,1),j=1,2)
 361         read (ithep_pdb,*,err=111,end=111) (polthet(j,i),j=0,3)
 362         read (ithep_pdb,*,err=111,end=111) (gthet(j,i),j=1,3)
 363         read (ithep_pdb,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
 364         sigc0(i)=sigc0(i)**2
 365       enddo
 366       do i=1,ntyp
 367       athet(1,i,1,-1)=athet(1,i,1,1)
 368       athet(2,i,1,-1)=athet(2,i,1,1)
 369       bthet(1,i,1,-1)=-bthet(1,i,1,1)
 370       bthet(2,i,1,-1)=-bthet(2,i,1,1)
 371       athet(1,i,-1,1)=-athet(1,i,1,1)
 372       athet(2,i,-1,1)=-athet(2,i,1,1)
 373       bthet(1,i,-1,1)=bthet(1,i,1,1)
 374       bthet(2,i,-1,1)=bthet(2,i,1,1)
 375       enddo
 376       do i=-ntyp,-1
 377       a0thet(i)=a0thet(-i)
 378       athet(1,i,-1,-1)=athet(1,-i,1,1)
 379       athet(2,i,-1,-1)=-athet(2,-i,1,1)
 380       bthet(1,i,-1,-1)=bthet(1,-i,1,1)
 381       bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
 382       athet(1,i,-1,1)=athet(1,-i,1,1)
 383       athet(2,i,-1,1)=-athet(2,-i,1,1)
 384       bthet(1,i,-1,1)=-bthet(1,-i,1,1)
 385       bthet(2,i,-1,1)=bthet(2,-i,1,1)
 386       athet(1,i,1,-1)=-athet(1,-i,1,1)
 387       athet(2,i,1,-1)=athet(2,-i,1,1)
 388       bthet(1,i,1,-1)=bthet(1,-i,1,1)
 389       bthet(2,i,1,-1)=-bthet(2,-i,1,1)
 390       theta0(i)=theta0(-i)
 391       sig0(i)=sig0(-i)
 392       sigc0(i)=sigc0(-i)
 393        do j=0,3
 394         polthet(j,i)=polthet(j,-i)
 395        enddo
 396        do j=1,3
 397          gthet(j,i)=gthet(j,-i)
 398        enddo
 399       enddo
 400       write (2,*) "End reading THETA_PDB"
 401       close (ithep_pdb)
 402 #endif
 403       close(ithep)
 404 #ifdef CRYST_SC
 405 C
 406 C Read the parameters of the probability distribution/energy expression
 407 C of the side chains.
 408 C
 409       do i=1,ntyp
 410         read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
 411         if (i.eq.10) then
 412           dsc_inv(i)=0.0D0
 413         else
 414           dsc_inv(i)=1.0D0/dsc(i)
 415         endif
 416         if (i.ne.10) then
 417         do j=1,nlob(i)
 418           do k=1,3
 419             do l=1,3
 420               blower(l,k,j)=0.0D0
 421             enddo
 422           enddo
 423         enddo
 424         bsc(1,i)=0.0D0
 425         read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),
 426      &    ((blower(k,l,1),l=1,k),k=1,3)
 427         do j=2,nlob(i)
 428           read (irotam,*,end=112,err=112) bsc(j,i)
 429           read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),
 430      &                                 ((blower(k,l,j),l=1,k),k=1,3)
 431         enddo
 432         do j=1,nlob(i)
 433           do k=1,3
 434             do l=1,k
 435               akl=0.0D0
 436               do m=1,3
 437                 akl=akl+blower(k,m,j)*blower(l,m,j)
 438               enddo
 439               gaussc(k,l,j,i)=akl
 440               gaussc(l,k,j,i)=akl
 441             enddo
 442           enddo
 443         enddo
 444         endif
 445       enddo
 446       close (irotam)
 447       if (lprint) then
 448         write (iout,'(/a)') 'Parameters of side-chain local geometry'
 449         do i=1,ntyp
 450           nlobi=nlob(i)
 451           if (nlobi.gt.0) then
 452             if (LaTeX) then
 453               write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
 454      &         ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
 455                write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
 456      &                             'log h',(bsc(j,i),j=1,nlobi)
 457                write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
 458      &        'x',((censc(k,j,i),k=1,3),j=1,nlobi)
 459               do k=1,3
 460                 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
 461      &                 ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
 462               enddo
 463             else
 464               write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
 465               write (iout,'(a,f10.4,4(16x,f10.4))')
 466      &                             'Center  ',(bsc(j,i),j=1,nlobi)
 467               write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),
 468      &           j=1,nlobi)
 469               write (iout,'(a)')
 470             endif
 471           endif
 472         enddo
 473       endif
 474 #else
 475 C
 476 C Read scrot parameters for potentials determined from all-atom AM1 calculations
 477 C added by Urszula Kozlowska 07/11/2007
 478 C
 479       do i=1,ntyp
 480         read (irotam,*,end=112,err=112)
 481        if (i.eq.10) then
 482          read (irotam,*,end=112,err=112)
 483        else
 484          do j=1,65
 485            read(irotam,*,end=112,err=112) sc_parmin(j,i)
 486          enddo
 487        endif
 488       enddo
 489 C
 490 C Read the parameters of the probability distribution/energy expression
 491 C of the side chains.
 492 C
 493       write (2,*) "Start reading ROTAM_PDB"
 494
 495       do i=1,ntyp
 496         read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
 497         if (i.eq.10) then
 498           dsc_inv(i)=0.0D0
 499         else
 500           dsc_inv(i)=1.0D0/dsc(i)
 501         endif
 502         if (i.ne.10) then
 503         do j=1,nlob(i)
 504           do k=1,3
 505             do l=1,3
 506               blower(l,k,j)=0.0D0
 507             enddo
 508           enddo
 509         enddo
 510         bsc(1,i)=0.0D0
 511         read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),
 512      &    ((blower(k,l,1),l=1,k),k=1,3)
 513         do j=2,nlob(i)
 514           read (irotam_pdb,*,end=112,err=112) bsc(j,i)
 515           read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),
 516      &                                 ((blower(k,l,j),l=1,k),k=1,3)
 517         enddo
 518         do j=1,nlob(i)
 519           do k=1,3
 520             do l=1,k
 521               akl=0.0D0
 522               do m=1,3
 523                 akl=akl+blower(k,m,j)*blower(l,m,j)
 524               enddo
 525               gaussc(k,l,j,i)=akl
 526               gaussc(l,k,j,i)=akl
 527             enddo
 528           enddo
 529         enddo
 530         endif
 531       enddo
 532       close (irotam_pdb)
 533       write (2,*) "Ending reading ROTAM_PDB"
 534 #endif
 535       close(irotam)
 536
 537 #ifdef CRYST_TOR
 538 C
 539 C Read torsional parameters in old format
 540 C
 541       read (itorp,*,end=113,err=113) ntortyp,nterm_old
 542       if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
 543       read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
 544       do i=1,ntortyp
 545         do j=1,ntortyp
 546           read (itorp,'(a)')
 547           do k=1,nterm_old
 548             read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
 549           enddo
 550         enddo
 551       enddo
 552       close (itorp)
 553       if (lprint) then
 554         write (iout,'(/a/)') 'Torsional constants:'
 555         do i=1,ntortyp
 556           do j=1,ntortyp
 557             write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
 558             write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
 559           enddo
 560         enddo
 561       endif
 562 #else
 563 C
 564 C Read torsional parameters
 565 C
 566       read (itorp,*,end=113,err=113) ntortyp
 567       read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
 568 c      write (iout,*) 'ntortyp',ntortyp
 569       do i=1,ntortyp
 570         do j=1,ntortyp
 571           read (itorp,*,end=113,err=113) nterm(i,j),nlor(i,j)
 572           v0ij=0.0d0
 573           si=-1.0d0
 574           do k=1,nterm(i,j)
 575             read (itorp,*,end=113,err=113) kk,v1(k,i,j),v2(k,i,j)
 576             v0ij=v0ij+si*v1(k,i,j)
 577             si=-si
 578           enddo
 579           do k=1,nlor(i,j)
 580             read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),
 581      &        vlor2(k,i,j),vlor3(k,i,j)
 582             v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
 583           enddo
 584           v0(i,j)=v0ij
 585         enddo
 586       enddo
 587       close (itorp)
 588       if (lprint) then
 589         write (iout,'(/a/)') 'Torsional constants:'
 590         do i=1,ntortyp
 591           do j=1,ntortyp
 592             write (iout,*) 'ityp',i,' jtyp',j
 593             write (iout,*) 'Fourier constants'
 594             do k=1,nterm(i,j)
 595               write (iout,'(2(1pe15.5))') v1(k,i,j),v2(k,i,j)
 596             enddo
 597             write (iout,*) 'Lorenz constants'
 598             do k=1,nlor(i,j)
 599               write (iout,'(3(1pe15.5))')
 600      &         vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
 601             enddo
 602           enddo
 603         enddo
 604       endif
 605 C
 606 C 6/23/01 Read parameters for double torsionals
 607 C
 608       do i=1,ntortyp
 609         do j=1,ntortyp
 610           do k=1,ntortyp
 611             read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
 612 c              write (iout,*) "OK onelett",
 613 c     &         i,j,k,t1,t2,t3
 614
 615             if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
 616      &        .or. t3.ne.toronelet(k)) then
 617               write (iout,*) "Error in double torsional parameter file",
 618      &         i,j,k,t1,t2,t3
 619 #ifdef MPI
 620               call MPI_Finalize(Ierror)
 621 #endif
 622                stop "Error in double torsional parameter file"
 623             endif
 624             read (itordp,*,end=114,err=114) ntermd_1(i,j,k),
 625      &         ntermd_2(i,j,k)
 626             read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k),l=1,
 627      &         ntermd_1(i,j,k))
 628             read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k),l=1,
 629      &         ntermd_1(i,j,k))
 630             read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k),l=1,
 631      &         ntermd_1(i,j,k))
 632             read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k),l=1,
 633      &         ntermd_1(i,j,k))
 634             read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k),
 635      &         v2c(m,l,i,j,k),v2s(l,m,i,j,k),v2s(m,l,i,j,k),
 636      &         m=1,l-1),l=1,ntermd_2(i,j,k))
 637           enddo
 638         enddo
 639       enddo
 640       if (lprint) then
 641       write (iout,*)
 642       write (iout,*) 'Constants for double torsionals'
 643       do iblock=1,2
 644       do i=0,ntortyp-1
 645         do j=-ntortyp+1,ntortyp-1
 646           do k=-ntortyp+1,ntortyp-1
 647             write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
 648      &        ' nsingle',ntermd_1(i,j,k),' ndouble',ntermd_2(i,j,k)
 649             write (iout,*)
 650             write (iout,*) 'Single angles:'
 651             do l=1,ntermd_1(i,j,k)
 652               write (iout,'(i5,2f10.5,5x,2f10.5)') l,
 653      &           v1c(1,l,i,j,k),v1s(1,l,i,j,k),
 654      &           v1c(2,l,i,j,k),v1s(2,l,i,j,k)
 655             enddo
 656             write (iout,*)
 657             write (iout,*) 'Pairs of angles:'
 658             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k))
 659             do l=1,ntermd_2(i,j,k)
 660               write (iout,'(i5,20f10.5)')
 661      &         l,(v2c(l,m,i,j,k),m=1,ntermd_2(i,j,k))
 662             enddo
 663             write (iout,*)
 664             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k))
 665             do l=1,ntermd_2(i,j,k)
 666               write (iout,'(i5,20f10.5)')
 667      &         l,(v2s(l,m,i,j,k),m=1,ntermd_2(i,j,k))
 668             enddo
 669             write (iout,*)
 670           enddo
 671         enddo
 672       enddo
 673       endif
 674 #endif
 675 C Read of Side-chain backbone correlation parameters
 676 C Modified 11 May 2012 by Adasko
 677 CCC
 678 C
 679       read (isccor,*,end=1113,err=1113) nsccortyp
 680 #ifdef SCCORPDB
 681       read (isccor,*,end=1113,err=1113) (isccortyp(i),i=1,ntyp)
 682       do i=-ntyp,-1
 683         isccortyp(i)=-isccortyp(-i)
 684       enddo
 685       iscprol=isccortyp(20)
 686 c      write (iout,*) 'ntortyp',ntortyp
 687       maxinter=3
 688 cc maxinter is maximum interaction sites
 689       do l=1,maxinter
 690       do i=1,nsccortyp
 691         do j=1,nsccortyp
 692           read (isccor,*,end=1113,err=1113) nterm_sccor(i,j),
 693      &             nlor_sccor(i,j)
 694           v0ijsccor=0.0d0
 695           v0ijsccor1=0.0d0
 696           v0ijsccor2=0.0d0
 697           v0ijsccor3=0.0d0
 698           si=-1.0d0
 699           nterm_sccor(-i,j)=nterm_sccor(i,j)
 700           nterm_sccor(-i,-j)=nterm_sccor(i,j)
 701           nterm_sccor(i,-j)=nterm_sccor(i,j)
 702           do k=1,nterm_sccor(i,j)
 703             read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
 704      &    ,v2sccor(k,l,i,j)
 705             if (j.eq.iscprol) then
 706               if (i.eq.isccortyp(10)) then
 707               v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
 708               v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
 709               else
 710              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
 711      &                        +v2sccor(k,l,i,j)*dsqrt(0.75d0)
 712              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
 713      &                        +v1sccor(k,l,i,j)*dsqrt(0.75d0)
 714              v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
 715              v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
 716              v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
 717              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
 718              endif
 719             else
 720               if (i.eq.isccortyp(10)) then
 721               v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
 722               v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
 723               else
 724                 if (j.eq.isccortyp(10)) then
 725               v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
 726               v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
 727                 else
 728              v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
 729              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
 730              v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
 731              v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
 732              v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
 733              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
 734             endif
 735              endif
 736              endif
 737             read (isccor,*,end=1113,err=1113) kk,v1sccor(k,l,i,j)
 738      &    ,v2sccor(k,l,i,j)
 739             v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
 740             v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
 741             v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
 742             v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
 743             si=-si
 744           enddo
 745           do k=1,nlor_sccor(i,j)
 746             read (isccor,*,end=1113,err=1113) kk,vlor1sccor(k,i,j),
 747      &        vlor2sccor(k,i,j),vlor3sccor(k,i,j)
 748             v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
 749      &(1+vlor3sccor(k,i,j)**2)
 750           enddo
 751           v0sccor(l,i,j)=v0ijsccor
 752           v0sccor(l,-i,j)=v0ijsccor1
 753           v0sccor(l,i,-j)=v0ijsccor2
 754           v0sccor(l,-i,-j)=v0ijsccor3
 755         enddo
 756       enddo
 757       enddo
 758       close (isccor)
 759 #else
 760       read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
 761 c      write (iout,*) 'ntortyp',ntortyp
 762       maxinter=3
 763 cc maxinter is maximum interaction sites
 764       do l=1,maxinter
 765       do i=1,nsccortyp
 766         do j=1,nsccortyp
 767           read (isccor,*,end=119,err=119)
 768      & nterm_sccor(i,j),nlor_sccor(i,j)
 769           v0ijsccor=0.0d0
 770           si=-1.0d0
 771
 772           do k=1,nterm_sccor(i,j)
 773             read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
 774      &    ,v2sccor(k,l,i,j)
 775             v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
 776             si=-si
 777           enddo
 778           do k=1,nlor_sccor(i,j)
 779             read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
 780      &        vlor2sccor(k,i,j),vlor3sccor(k,i,j)
 781             v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
 782      &(1+vlor3sccor(k,i,j)**2)
 783           enddo
 784           v0sccor(i,j,iblock)=v0ijsccor
 785         enddo
 786       enddo
 787       enddo
 788       close (isccor)
 789
 790 #endif
 791       if (lprint) then
 792         write (iout,'(/a/)') 'Torsional constants:'
 793         do i=1,nsccortyp
 794           do j=1,nsccortyp
 795             write (iout,*) 'ityp',i,' jtyp',j
 796             write (iout,*) 'Fourier constants'
 797             do k=1,nterm_sccor(i,j)
 798       write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
 799             enddo
 800             write (iout,*) 'Lorenz constants'
 801             do k=1,nlor_sccor(i,j)
 802               write (iout,'(3(1pe15.5))')
 803      &         vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
 804             enddo
 805           enddo
 806         enddo
 807       endif
 808 C
 809 C
 810 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
 811 C         interaction energy of the Gly, Ala, and Pro prototypes.
 812 C
 813       if (lprint) then
 814         write (iout,*)
 815         write (iout,*) "Coefficients of the cumulants"
 816       endif
 817       read (ifourier,*) nloctyp
 818       do i=1,nloctyp
 819         read (ifourier,*,end=115,err=115)
 820         read (ifourier,*,end=115,err=115) (b(ii),ii=1,13)
 821         if (lprint) then
 822         write (iout,*) 'Type',i
 823         write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
 824         endif
 825         B1(1,i)  = b(3)
 826         B1(2,i)  = b(5)
 827 c        b1(1,i)=0.0d0
 828 c        b1(2,i)=0.0d0
 829         B1tilde(1,i) = b(3)
 830         B1tilde(2,i) =-b(5)
 831         B1tilde(1,-i) =-b(3)
 832         B1tilde(2,-i) =b(5)
 833 c        b1tilde(1,i)=0.0d0
 834 c        b1tilde(2,i)=0.0d0
 835         B2(1,i)  = b(2)
 836         B2(2,i)  = b(4)
 837 c        b2(1,i)=0.0d0
 838 c        b2(2,i)=0.0d0
 839         CC(1,1,i)= b(7)
 840         CC(2,2,i)=-b(7)
 841         CC(2,1,i)= b(9)
 842         CC(1,2,i)= b(9)
 843 c        CC(1,1,i)=0.0d0
 844 c        CC(2,2,i)=0.0d0
 845 c        CC(2,1,i)=0.0d0
 846 c        CC(1,2,i)=0.0d0
 847         Ctilde(1,1,i)=b(7)
 848         Ctilde(1,2,i)=b(9)
 849         Ctilde(2,1,i)=-b(9)
 850         Ctilde(2,2,i)=b(7)
 851 c        Ctilde(1,1,i)=0.0d0
 852 c        Ctilde(1,2,i)=0.0d0
 853 c        Ctilde(2,1,i)=0.0d0
 854 c        Ctilde(2,2,i)=0.0d0
 855         DD(1,1,i)= b(6)
 856         DD(2,2,i)=-b(6)
 857         DD(2,1,i)= b(8)
 858         DD(1,2,i)= b(8)
 859 c        DD(1,1,i)=0.0d0
 860 c        DD(2,2,i)=0.0d0
 861 c        DD(2,1,i)=0.0d0
 862 c        DD(1,2,i)=0.0d0
 863         Dtilde(1,1,i)=b(6)
 864         Dtilde(1,2,i)=b(8)
 865         Dtilde(2,1,i)=-b(8)
 866         Dtilde(2,2,i)=b(6)
 867 c        Dtilde(1,1,i)=0.0d0
 868 c        Dtilde(1,2,i)=0.0d0
 869 c        Dtilde(2,1,i)=0.0d0
 870 c        Dtilde(2,2,i)=0.0d0
 871         EE(1,1,i)= b(10)+b(11)
 872         EE(2,2,i)=-b(10)+b(11)
 873         EE(2,1,i)= b(12)-b(13)
 874         EE(1,2,i)= b(12)+b(13)
 875 c        ee(1,1,i)=1.0d0
 876 c        ee(2,2,i)=1.0d0
 877 c        ee(2,1,i)=0.0d0
 878 c        ee(1,2,i)=0.0d0
 879 c        ee(2,1,i)=ee(1,2,i)
 880       enddo
 881       if (lprint) then
 882       do i=1,nloctyp
 883         write (iout,*) 'Type',i
 884         write (iout,*) 'B1'
 885         write(iout,*) B1(1,i),B1(2,i)
 886         write (iout,*) 'B2'
 887         write(iout,*) B2(1,i),B2(2,i)
 888         write (iout,*) 'CC'
 889         do j=1,2
 890           write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
 891         enddo
 892         write(iout,*) 'DD'
 893         do j=1,2
 894           write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
 895         enddo
 896         write(iout,*) 'EE'
 897         do j=1,2
 898           write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
 899         enddo
 900       enddo
 901       endif
 902 C
 903 C Read electrostatic-interaction parameters
 904 C
 905       if (lprint) then
 906         write (iout,*)
 907         write (iout,'(/a)') 'Electrostatic interaction constants:'
 908         write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
 909      &            'IT','JT','APP','BPP','AEL6','AEL3'
 910       endif
 911       read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
 912       read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
 913       read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
 914       read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
 915       close (ielep)
 916       do i=1,2
 917         do j=1,2
 918         rri=rpp(i,j)**6
 919         app (i,j)=epp(i,j)*rri*rri
 920         bpp (i,j)=-2.0D0*epp(i,j)*rri
 921         ael6(i,j)=elpp6(i,j)*4.2D0**6
 922         ael3(i,j)=elpp3(i,j)*4.2D0**3
 923         if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
 924      &                    ael6(i,j),ael3(i,j)
 925         enddo
 926       enddo
 927 C
 928 C Read side-chain interaction parameters.
 929 C
 930       read (isidep,*,end=117,err=117) ipot,expon
 931       if (ipot.lt.1 .or. ipot.gt.5) then
 932         write (iout,'(2a)') 'Error while reading SC interaction',
 933      &               'potential file - unknown potential type.'
 934 #ifdef MPI
 935         call MPI_Finalize(Ierror)
 936 #endif
 937         stop
 938       endif
 939       expon2=expon/2
 940       if(me.eq.king)
 941      & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
 942      & ', exponents are ',expon,2*expon
 943       goto (10,20,30,30,40) ipot
 944 C----------------------- LJ potential ---------------------------------
 945    10 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
 946      &   (sigma0(i),i=1,ntyp)
 947       if (lprint) then
 948         write (iout,'(/a/)') 'Parameters of the LJ potential:'
 949         write (iout,'(a/)') 'The epsilon array:'
 950         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 951         write (iout,'(/a)') 'One-body parameters:'
 952         write (iout,'(a,4x,a)') 'residue','sigma'
 953         write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
 954       endif
 955       goto 50
 956 C----------------------- LJK potential --------------------------------
 957    20 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
 958      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
 959       if (lprint) then
 960         write (iout,'(/a/)') 'Parameters of the LJK potential:'
 961         write (iout,'(a/)') 'The epsilon array:'
 962         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 963         write (iout,'(/a)') 'One-body parameters:'
 964         write (iout,'(a,4x,2a)') 'residue','   sigma  ','    r0    '
 965         write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
 966      &        i=1,ntyp)
 967       endif
 968       goto 50
 969 C---------------------- GB or BP potential -----------------------------
 970    30 do i=1,ntyp
 971         read (isidep,*,end=116,err=116)(eps(i,j),j=i,ntyp)
 972        enddo
 973        read (isidep,*,end=116,err=116)(sigma0(i),i=1,ntyp)
 974        read (isidep,*,end=116,err=116)(sigii(i),i=1,ntyp)
 975        read (isidep,*,end=116,err=116)(chip(i),i=1,ntyp)
 976        read (isidep,*,end=116,err=116)(alp(i),i=1,ntyp)
 977
 978 c   30 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
 979 c     &  (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip(i),i=1,ntyp),
 980 c     &  (alp(i),i=1,ntyp)
 981 C For the GB potential convert sigma'**2 into chi'
 982       if (ipot.eq.4) then
 983         do i=1,ntyp
 984           chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
 985         enddo
 986       endif
 987       if (lprint) then
 988         write (iout,'(/a/)') 'Parameters of the BP potential:'
 989         write (iout,'(a/)') 'The epsilon array:'
 990         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 991         write (iout,'(/a)') 'One-body parameters:'
 992         write (iout,'(a,4x,4a)') 'residue','   sigma  ','s||/s_|_^2',
 993      &       '    chip  ','    alph  '
 994         write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
 995      &                     chip(i),alp(i),i=1,ntyp)
 996       endif
 997       goto 50
 998 C--------------------- GBV potential -----------------------------------
 999    40 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
1000      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1001      &  (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1002       if (lprint) then
1003         write (iout,'(/a/)') 'Parameters of the GBV potential:'
1004         write (iout,'(a/)') 'The epsilon array:'
1005         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1006         write (iout,'(/a)') 'One-body parameters:'
1007         write (iout,'(a,4x,5a)') 'residue','   sigma  ','    r0    ',
1008      &      's||/s_|_^2','    chip  ','    alph  '
1009         write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1010      &           sigii(i),chip(i),alp(i),i=1,ntyp)
1011       endif
1012    50 continue
1013       close (isidep)
1014 C-----------------------------------------------------------------------
1015 C Calculate the "working" parameters of SC interactions.
1016       do i=2,ntyp
1017         do j=1,i-1
1018           eps(i,j)=eps(j,i)
1019         enddo
1020       enddo
1021       do i=1,ntyp
1022         do j=i,ntyp
1023           sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1024           sigma(j,i)=sigma(i,j)
1025           rs0(i,j)=dwa16*sigma(i,j)
1026           rs0(j,i)=rs0(i,j)
1027         enddo
1028       enddo
1029       if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1030      & 'Working parameters of the SC interactions:',
1031      & '     a    ','     b    ','   augm   ','  sigma ','   r0   ',
1032      & '  chi1   ','   chi2   '
1033       do i=1,ntyp
1034         do j=i,ntyp
1035           epsij=eps(i,j)
1036           if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1037             rrij=sigma(i,j)
1038           else
1039             rrij=rr0(i)+rr0(j)
1040           endif
1041           r0(i,j)=rrij
1042           r0(j,i)=rrij
1043           rrij=rrij**expon
1044           epsij=eps(i,j)
1045           sigeps=dsign(1.0D0,epsij)
1046           epsij=dabs(epsij)
1047           aa(i,j)=epsij*rrij*rrij
1048           bb(i,j)=-sigeps*epsij*rrij
1049           aa(j,i)=aa(i,j)
1050           bb(j,i)=bb(i,j)
1051           if (ipot.gt.2) then
1052             sigt1sq=sigma0(i)**2
1053             sigt2sq=sigma0(j)**2
1054             sigii1=sigii(i)
1055             sigii2=sigii(j)
1056             ratsig1=sigt2sq/sigt1sq
1057             ratsig2=1.0D0/ratsig1
1058             chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1059             if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1060             rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1061           else
1062             rsum_max=sigma(i,j)
1063           endif
1064 c         if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1065             sigmaii(i,j)=rsum_max
1066             sigmaii(j,i)=rsum_max
1067 c         else
1068 c           sigmaii(i,j)=r0(i,j)
1069 c           sigmaii(j,i)=r0(i,j)
1070 c         endif
1071 cd        write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1072           if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1073             r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1074             augm(i,j)=epsij*r_augm**(2*expon)
1075 c           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1076             augm(j,i)=augm(i,j)
1077           else
1078             augm(i,j)=0.0D0
1079             augm(j,i)=0.0D0
1080           endif
1081           if (lprint) then
1082             write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1083      &      restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1084      &      sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1085           endif
1086         enddo
1087       enddo
1088 #ifdef OLDSCP
1089 C
1090 C Define the SC-p interaction constants (hard-coded; old style)
1091 C
1092       do i=1,ntyp
1093 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1094 C helix formation)
1095 c       aad(i,1)=0.3D0*4.0D0**12
1096 C Following line for constants currently implemented
1097 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1098         aad(i,1)=1.5D0*4.0D0**12
1099 c       aad(i,1)=0.17D0*5.6D0**12
1100         aad(i,2)=aad(i,1)
1101 C "Soft" SC-p repulsion
1102         bad(i,1)=0.0D0
1103 C Following line for constants currently implemented
1104 c       aad(i,1)=0.3D0*4.0D0**6
1105 C "Hard" SC-p repulsion
1106         bad(i,1)=3.0D0*4.0D0**6
1107 c       bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1108         bad(i,2)=bad(i,1)
1109 c       aad(i,1)=0.0D0
1110 c       aad(i,2)=0.0D0
1111 c       bad(i,1)=1228.8D0
1112 c       bad(i,2)=1228.8D0
1113       enddo
1114 #else
1115 C
1116 C 8/9/01 Read the SC-p interaction constants from file
1117 C
1118       do i=1,ntyp
1119         read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1120       enddo
1121       do i=1,ntyp
1122         aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1123         aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1124         bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1125         bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1126       enddo
1127
1128       if (lprint) then
1129         write (iout,*) "Parameters of SC-p interactions:"
1130         do i=1,ntyp
1131           write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1132      &     eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1133         enddo
1134       endif
1135 #endif
1136 C
1137 C Define the constants of the disulfide bridge
1138 C
1139       ebr=-5.50D0
1140 c
1141 c Old arbitrary potential - commented out.
1142 c
1143 c      dbr= 4.20D0
1144 c      fbr= 3.30D0
1145 c
1146 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1147 c energy surface of diethyl disulfide.
1148 c A. Liwo and U. Kozlowska, 11/24/03
1149 c
1150       D0CM = 3.78d0
1151       AKCM = 15.1d0
1152       AKTH = 11.0d0
1153       AKCT = 12.0d0
1154       V1SS =-1.08d0
1155       V2SS = 7.61d0
1156       V3SS = 13.7d0
1157 c      akcm=0.0d0
1158 c      akth=0.0d0
1159 c      akct=0.0d0
1160 c      v1ss=0.0d0
1161 c      v2ss=0.0d0
1162 c      v3ss=0.0d0
1163
1164       if(me.eq.king) then
1165       write (iout,'(/a)') "Disulfide bridge parameters:"
1166       write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1167       write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1168       write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1169       write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1170      &  ' v3ss:',v3ss
1171       endif
1172       return
1173   111 write (iout,*) "Error reading bending energy parameters."
1174       goto 999
1175   112 write (iout,*) "Error reading rotamer energy parameters."
1176       goto 999
1177   113 write (iout,*) "Error reading torsional energy parameters."
1178       goto 999
1179  1113 write (iout,*)
1180      &  "Error reading side-chain torsional energy parameters."
1181       goto 999
1182   114 write (iout,*) "Error reading double torsional energy parameters."
1183       goto 999
1184   115 write (iout,*)
1185      &  "Error reading cumulant (multibody energy) parameters."
1186       goto 999
1187   116 write (iout,*) "Error reading electrostatic energy parameters."
1188       goto 999
1189   117 write (iout,*) "Error reading side chain interaction parameters."
1190       goto 999
1191   118 write (iout,*) "Error reading SCp interaction parameters."
1192       goto 999
1193   119 write (iout,*) "Error reading SCCOR parameters"
1194   999 continue
1195 #ifdef MPI
1196       call MPI_Finalize(Ierror)
1197 #endif
1198       stop
1199       return
1200       end
1201
1202
1203       subroutine getenv_loc(var, val)
1204       character(*) var, val
1205
1206 #ifdef WINIFL
1207       character(2000) line
1208       external ilen
1209
1210       open (196,file='env',status='old',readonly,shared)
1211       iread=0
1212 c      write(*,*)'looking for ',var
1213 10    read(196,*,err=11,end=11)line
1214       iread=index(line,var)
1215 c      write(*,*)iread,' ',var,' ',line
1216       if (iread.eq.0) go to 10
1217 c      write(*,*)'---> ',line
1218 11    continue
1219       if(iread.eq.0) then
1220 c       write(*,*)'CHUJ'
1221        val=''
1222       else
1223        iread=iread+ilen(var)+1
1224        read (line(iread:),*,err=12,end=12) val
1225 c       write(*,*)'OK: ',var,' = ',val
1226       endif
1227       close(196)
1228       return
1229 12    val=''
1230       close(196)
1231 #elif (defined CRAY)
1232       integer lennam,lenval,ierror
1233 c
1234 c        getenv using a POSIX call, useful on the T3D
1235 c        Sept 1996, comment out error check on advice of H. Pritchard
1236 c
1237       lennam = len(var)
1238       if(lennam.le.0) stop '--error calling getenv--'
1239       call pxfgetenv(var,lennam,val,lenval,ierror)
1240 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1241 #else
1242       call getenv(var,val)
1243 #endif
1244
1245       return
1246       end