source/unres/src_MD/parmread.F

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