source/unres/src_MD-M/parmread.F

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