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