source/unres/src_MD-M/parmread.F

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