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           si=-1.0d0
 624           nterm_sccor(-i,j)=nterm_sccor(i,j)
 625           nterm_sccor(-i,-j)=nterm_sccor(i,j)
 626           nterm_sccor(i,-j)=nterm_sccor(i,j)
 627           do k=1,nterm_sccor(i,j)
 628             read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
 629      &    ,v2sccor(k,l,i,j)
 630             if (j.eq.iscprol) then
 631              if (i.eq.isccortyp(10)) then
 632              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
 633              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
 634              else
 635              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
 636      &                        +v2sccor(k,l,i,j)*dsqrt(0.75d0)
 637              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
 638      &                        +v1sccor(k,l,i,j)*dsqrt(0.75d0)
 639              v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
 640              v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
 641              v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
 642              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
 643              endif
 644             else
 645              if (i.eq.isccortyp(10)) then
 646              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
 647              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
 648              else
 649                if (j.eq.isccortyp(10)) then
 650              v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
 651              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
 652                else
 653              v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
 654              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
 655              v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
 656              v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
 657              v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
 658              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
 659                 endif
 660                endif
 661             endif
 662             v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
 663             si=-si
 664           enddo
 665           do k=1,nlor_sccor(i,j)
 666             read (isccor,*,end=113,err=113) kk,vlor1sccor(k,i,j),
 667      &        vlor2sccor(k,i,j),vlor3sccor(k,i,j)
 668             v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
 669      &(1+vlor3sccor(k,i,j)**2)
 670           enddo
 671           v0sccor(i,j)=v0ijsccor
 672         enddo
 673       enddo
 674       enddo
 675       close (isccor)
 676 #else
 677       read (isccor,*,end=113,err=113) (isccortyp(i),i=1,ntyp)
 678 c      write (iout,*) 'ntortyp',ntortyp
 679       maxinter=3
 680 cc maxinter is maximum interaction sites
 681       do l=1,maxinter
 682       do i=1,nsccortyp
 683         do j=1,nsccortyp
 684           read (isccor,*,end=113,err=113)
 685      & nterm_sccor(i,j),nlor_sccor(i,j)
 686           v0ijsccor=0.0d0
 687           si=-1.0d0
 688
 689           do k=1,nterm_sccor(i,j)
 690             read (isccor,*,end=113,err=113) kk,v1sccor(k,l,i,j)
 691      &    ,v2sccor(k,l,i,j)
 692             v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
 693             si=-si
 694           enddo
 695           do k=1,nlor_sccor(i,j)
 696             read (isccor,*,end=113,err=113) kk,vlor1sccor(k,i,j),
 697      &        vlor2sccor(k,i,j),vlor3sccor(k,i,j)
 698             v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
 699      &(1+vlor3sccor(k,i,j)**2)
 700           enddo
 701           v0sccor(i,j)=v0ijsccor
 702         enddo
 703       enddo
 704       enddo
 705       close (isccor)
 706
 707 #endif
 708       if (lprint) then
 709         write (iout,'(/a/)') 'Torsional constants:'
 710         do i=1,nsccortyp
 711           do j=1,nsccortyp
 712             write (iout,*) 'ityp',i,' jtyp',j
 713             write (iout,*) 'Fourier constants'
 714             do k=1,nterm_sccor(i,j)
 715       write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
 716             enddo
 717             write (iout,*) 'Lorenz constants'
 718             do k=1,nlor_sccor(i,j)
 719               write (iout,'(3(1pe15.5))')
 720      &         vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
 721             enddo
 722           enddo
 723         enddo
 724       endif
 725
 726 C
 727 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
 728 C         interaction energy of the Gly, Ala, and Pro prototypes.
 729 C
 730       if (lprint) then
 731         write (iout,*)
 732         write (iout,*) "Coefficients of the cumulants"
 733       endif
 734       read (ifourier,*) nloctyp
 735       do i=0,nloctyp-1
 736         read (ifourier,*,end=115,err=115)
 737         read (ifourier,*,end=115,err=115) (b(ii),ii=1,13)
 738         if (lprint) then
 739         write (iout,*) 'Type',i
 740         write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
 741         endif
 742         B1(1,i)  = b(3)
 743         B1(2,i)  = b(5)
 744         B1(1,-i) = b(3)
 745         B1(2,-i) = -b(5)
 746 c        b1(1,i)=0.0d0
 747 c        b1(2,i)=0.0d0
 748         B1tilde(1,i) = b(3)
 749         B1tilde(2,i) =-b(5)
 750         B1tilde(1,-i) =-b(3)
 751         B1tilde(2,-i) =b(5)
 752 c        b1tilde(1,i)=0.0d0
 753 c        b1tilde(2,i)=0.0d0
 754         B2(1,i)  = b(2)
 755         B2(2,i)  = b(4)
 756         B2(1,-i)  =b(2)
 757         B2(2,-i)  =-b(4)
 758
 759 c        b2(1,i)=0.0d0
 760 c        b2(2,i)=0.0d0
 761         CC(1,1,i)= b(7)
 762         CC(2,2,i)=-b(7)
 763         CC(2,1,i)= b(9)
 764         CC(1,2,i)= b(9)
 765         CC(1,1,-i)= b(7)
 766         CC(2,2,-i)=-b(7)
 767         CC(2,1,-i)=-b(9)
 768         CC(1,2,-i)=-b(9)
 769 c        CC(1,1,i)=0.0d0
 770 c        CC(2,2,i)=0.0d0
 771 c        CC(2,1,i)=0.0d0
 772 c        CC(1,2,i)=0.0d0
 773         Ctilde(1,1,i)=b(7)
 774         Ctilde(1,2,i)=b(9)
 775         Ctilde(2,1,i)=-b(9)
 776         Ctilde(2,2,i)=b(7)
 777         Ctilde(1,1,-i)=b(7)
 778         Ctilde(1,2,-i)=-b(9)
 779         Ctilde(2,1,-i)=b(9)
 780         Ctilde(2,2,-i)=b(7)
 781
 782 c        Ctilde(1,1,i)=0.0d0
 783 c        Ctilde(1,2,i)=0.0d0
 784 c        Ctilde(2,1,i)=0.0d0
 785 c        Ctilde(2,2,i)=0.0d0
 786         DD(1,1,i)= b(6)
 787         DD(2,2,i)=-b(6)
 788         DD(2,1,i)= b(8)
 789         DD(1,2,i)= b(8)
 790         DD(1,1,-i)= b(6)
 791         DD(2,2,-i)=-b(6)
 792         DD(2,1,-i)=-b(8)
 793         DD(1,2,-i)=-b(8)
 794 c        DD(1,1,i)=0.0d0
 795 c        DD(2,2,i)=0.0d0
 796 c        DD(2,1,i)=0.0d0
 797 c        DD(1,2,i)=0.0d0
 798         Dtilde(1,1,i)=b(6)
 799         Dtilde(1,2,i)=b(8)
 800         Dtilde(2,1,i)=-b(8)
 801         Dtilde(2,2,i)=b(6)
 802         Dtilde(1,1,-i)=b(6)
 803         Dtilde(1,2,-i)=-b(8)
 804         Dtilde(2,1,-i)=b(8)
 805         Dtilde(2,2,-i)=b(6)
 806
 807 c        Dtilde(1,1,i)=0.0d0
 808 c        Dtilde(1,2,i)=0.0d0
 809 c        Dtilde(2,1,i)=0.0d0
 810 c        Dtilde(2,2,i)=0.0d0
 811         EE(1,1,i)= b(10)+b(11)
 812         EE(2,2,i)=-b(10)+b(11)
 813         EE(2,1,i)= b(12)-b(13)
 814         EE(1,2,i)= b(12)+b(13)
 815         EE(1,1,-i)= b(10)+b(11)
 816         EE(2,2,-i)=-b(10)+b(11)
 817         EE(2,1,-i)=-b(12)+b(13)
 818         EE(1,2,-i)=-b(12)-b(13)
 819
 820 c        ee(1,1,i)=1.0d0
 821 c        ee(2,2,i)=1.0d0
 822 c        ee(2,1,i)=0.0d0
 823 c        ee(1,2,i)=0.0d0
 824 c        ee(2,1,i)=ee(1,2,i)
 825       enddo
 826       if (lprint) then
 827       do i=1,nloctyp
 828         write (iout,*) 'Type',i
 829         write (iout,*) 'B1'
 830         write(iout,*) B1(1,i),B1(2,i)
 831         write (iout,*) 'B2'
 832         write(iout,*) B2(1,i),B2(2,i)
 833         write (iout,*) 'CC'
 834         do j=1,2
 835           write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
 836         enddo
 837         write(iout,*) 'DD'
 838         do j=1,2
 839           write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
 840         enddo
 841         write(iout,*) 'EE'
 842         do j=1,2
 843           write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
 844         enddo
 845       enddo
 846       endif
 847
 848 C
 849 C Read electrostatic-interaction parameters
 850 C
 851       if (lprint) then
 852         write (iout,*)
 853         write (iout,'(/a)') 'Electrostatic interaction constants:'
 854         write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
 855      &            'IT','JT','APP','BPP','AEL6','AEL3'
 856       endif
 857       read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
 858       read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
 859       read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
 860       read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
 861       close (ielep)
 862       do i=1,2
 863         do j=1,2
 864         rri=rpp(i,j)**6
 865         app (i,j)=epp(i,j)*rri*rri
 866         bpp (i,j)=-2.0D0*epp(i,j)*rri
 867         ael6(i,j)=elpp6(i,j)*4.2D0**6
 868         ael3(i,j)=elpp3(i,j)*4.2D0**3
 869         if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
 870      &                    ael6(i,j),ael3(i,j)
 871         enddo
 872       enddo
 873 C
 874 C Read side-chain interaction parameters.
 875 C
 876       read (isidep,*,end=117,err=117) ipot,expon
 877       if (ipot.lt.1 .or. ipot.gt.5) then
 878         write (iout,'(2a)') 'Error while reading SC interaction',
 879      &               'potential file - unknown potential type.'
 880 #ifdef MPI
 881         call MPI_Finalize(Ierror)
 882 #endif
 883         stop
 884       endif
 885       expon2=expon/2
 886       if(me.eq.king)
 887      & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
 888      & ', exponents are ',expon,2*expon
 889       goto (10,20,30,30,40) ipot
 890 C----------------------- LJ potential ---------------------------------
 891    10 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
 892      &   (sigma0(i),i=1,ntyp)
 893       if (lprint) then
 894         write (iout,'(/a/)') 'Parameters of the LJ potential:'
 895         write (iout,'(a/)') 'The epsilon array:'
 896         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 897         write (iout,'(/a)') 'One-body parameters:'
 898         write (iout,'(a,4x,a)') 'residue','sigma'
 899         write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
 900       endif
 901       goto 50
 902 C----------------------- LJK potential --------------------------------
 903    20 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
 904      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
 905       if (lprint) then
 906         write (iout,'(/a/)') 'Parameters of the LJK potential:'
 907         write (iout,'(a/)') 'The epsilon array:'
 908         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 909         write (iout,'(/a)') 'One-body parameters:'
 910         write (iout,'(a,4x,2a)') 'residue','   sigma  ','    r0    '
 911         write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
 912      &        i=1,ntyp)
 913       endif
 914       goto 50
 915 C---------------------- GB or BP potential -----------------------------
 916    30 do i=1,ntyp
 917        read (isidep,*,end=116,err=116)(eps(i,j),j=i,ntyp)
 918       enddo
 919       read (isidep,*,end=116,err=116)(sigma0(i),i=1,ntyp)
 920       read (isidep,*,end=116,err=116)(sigii(i),i=1,ntyp)
 921       read (isidep,*,end=116,err=116)(chip(i),i=1,ntyp)
 922       read (isidep,*,end=116,err=116)(alp(i),i=1,ntyp)
 923 C For the GB potential convert sigma'**2 into chi'
 924       if (ipot.eq.4) then
 925         do i=1,ntyp
 926           chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
 927         enddo
 928       endif
 929       if (lprint) then
 930         write (iout,'(/a/)') 'Parameters of the BP potential:'
 931         write (iout,'(a/)') 'The epsilon array:'
 932         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 933         write (iout,'(/a)') 'One-body parameters:'
 934         write (iout,'(a,4x,4a)') 'residue','   sigma  ','s||/s_|_^2',
 935      &       '    chip  ','    alph  '
 936         write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
 937      &                     chip(i),alp(i),i=1,ntyp)
 938       endif
 939       goto 50
 940 C--------------------- GBV potential -----------------------------------
 941    40 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
 942      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
 943      &  (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
 944       if (lprint) then
 945         write (iout,'(/a/)') 'Parameters of the GBV potential:'
 946         write (iout,'(a/)') 'The epsilon array:'
 947         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 948         write (iout,'(/a)') 'One-body parameters:'
 949         write (iout,'(a,4x,5a)') 'residue','   sigma  ','    r0    ',
 950      &      's||/s_|_^2','    chip  ','    alph  '
 951         write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
 952      &           sigii(i),chip(i),alp(i),i=1,ntyp)
 953       endif
 954    50 continue
 955       close (isidep)
 956 C-----------------------------------------------------------------------
 957 C Calculate the "working" parameters of SC interactions.
 958       do i=2,ntyp
 959         do j=1,i-1
 960           eps(i,j)=eps(j,i)
 961         enddo
 962       enddo
 963       do i=1,ntyp
 964         do j=i,ntyp
 965           sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
 966           sigma(j,i)=sigma(i,j)
 967           rs0(i,j)=dwa16*sigma(i,j)
 968           rs0(j,i)=rs0(i,j)
 969         enddo
 970       enddo
 971       if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
 972      & 'Working parameters of the SC interactions:',
 973      & '     a    ','     b    ','   augm   ','  sigma ','   r0   ',
 974      & '  chi1   ','   chi2   '
 975       do i=1,ntyp
 976         do j=i,ntyp
 977           epsij=eps(i,j)
 978           if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
 979             rrij=sigma(i,j)
 980           else
 981             rrij=rr0(i)+rr0(j)
 982           endif
 983           r0(i,j)=rrij
 984           r0(j,i)=rrij
 985           rrij=rrij**expon
 986           epsij=eps(i,j)
 987           sigeps=dsign(1.0D0,epsij)
 988           epsij=dabs(epsij)
 989           aa(i,j)=epsij*rrij*rrij
 990           bb(i,j)=-sigeps*epsij*rrij
 991           aa(j,i)=aa(i,j)
 992           bb(j,i)=bb(i,j)
 993           if (ipot.gt.2) then
 994             sigt1sq=sigma0(i)**2
 995             sigt2sq=sigma0(j)**2
 996             sigii1=sigii(i)
 997             sigii2=sigii(j)
 998             ratsig1=sigt2sq/sigt1sq
 999             ratsig2=1.0D0/ratsig1
1000             chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1001             if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1002             rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1003           else
1004             rsum_max=sigma(i,j)
1005           endif
1006 c         if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1007             sigmaii(i,j)=rsum_max
1008             sigmaii(j,i)=rsum_max
1009 c         else
1010 c           sigmaii(i,j)=r0(i,j)
1011 c           sigmaii(j,i)=r0(i,j)
1012 c         endif
1013 cd        write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1014           if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1015             r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1016             augm(i,j)=epsij*r_augm**(2*expon)
1017 c           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1018             augm(j,i)=augm(i,j)
1019           else
1020             augm(i,j)=0.0D0
1021             augm(j,i)=0.0D0
1022           endif
1023           if (lprint) then
1024             write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1025      &      restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1026      &      sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1027           endif
1028         enddo
1029       enddo
1030 #ifdef OLDSCP
1031 C
1032 C Define the SC-p interaction constants (hard-coded; old style)
1033 C
1034       do i=1,ntyp
1035 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1036 C helix formation)
1037 c       aad(i,1)=0.3D0*4.0D0**12
1038 C Following line for constants currently implemented
1039 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1040         aad(i,1)=1.5D0*4.0D0**12
1041 c       aad(i,1)=0.17D0*5.6D0**12
1042         aad(i,2)=aad(i,1)
1043 C "Soft" SC-p repulsion
1044         bad(i,1)=0.0D0
1045 C Following line for constants currently implemented
1046 c       aad(i,1)=0.3D0*4.0D0**6
1047 C "Hard" SC-p repulsion
1048         bad(i,1)=3.0D0*4.0D0**6
1049 c       bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1050         bad(i,2)=bad(i,1)
1051 c       aad(i,1)=0.0D0
1052 c       aad(i,2)=0.0D0
1053 c       bad(i,1)=1228.8D0
1054 c       bad(i,2)=1228.8D0
1055       enddo
1056 #else
1057 C
1058 C 8/9/01 Read the SC-p interaction constants from file
1059 C
1060       do i=1,ntyp
1061         read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1062       enddo
1063       do i=1,ntyp
1064         aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1065         aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1066         bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1067         bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1068       enddo
1069
1070       if (lprint) then
1071         write (iout,*) "Parameters of SC-p interactions:"
1072         do i=1,ntyp
1073           write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1074      &     eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1075         enddo
1076       endif
1077 #endif
1078 C
1079 C Define the constants of the disulfide bridge
1080 C
1081       ebr=-5.50D0
1082 c
1083 c Old arbitrary potential - commented out.
1084 c
1085 c      dbr= 4.20D0
1086 c      fbr= 3.30D0
1087 c
1088 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1089 c energy surface of diethyl disulfide.
1090 c A. Liwo and U. Kozlowska, 11/24/03
1091 c
1092       D0CM = 3.78d0
1093       AKCM = 15.1d0
1094       AKTH = 11.0d0
1095       AKCT = 12.0d0
1096       V1SS =-1.08d0
1097       V2SS = 7.61d0
1098       V3SS = 13.7d0
1099 c      akcm=0.0d0
1100 c      akth=0.0d0
1101 c      akct=0.0d0
1102 c      v1ss=0.0d0
1103 c      v2ss=0.0d0
1104 c      v3ss=0.0d0
1105
1106       if(me.eq.king) then
1107       write (iout,'(/a)') "Disulfide bridge parameters:"
1108       write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1109       write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1110       write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1111       write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1112      &  ' v3ss:',v3ss
1113       endif
1114       return
1115   111 write (iout,*) "Error reading bending energy parameters."
1116       goto 999
1117   112 write (iout,*) "Error reading rotamer energy parameters."
1118       goto 999
1119   113 write (iout,*) "Error reading torsional energy parameters."
1120       goto 999
1121   114 write (iout,*) "Error reading double torsional energy parameters."
1122       goto 999
1123   115 write (iout,*)
1124      &  "Error reading cumulant (multibody energy) parameters."
1125       goto 999
1126   116 write (iout,*) "Error reading electrostatic energy parameters."
1127       goto 999
1128   117 write (iout,*) "Error reading side chain interaction parameters."
1129       goto 999
1130   118 write (iout,*) "Error reading SCp interaction parameters."
1131       goto 999
1132   119 write (iout,*) "Error reading SCCOR parameters"
1133   999 continue
1134 #ifdef MPI
1135       call MPI_Finalize(Ierror)
1136 #endif
1137       stop
1138       return
1139       end
1140
1141
1142       subroutine getenv_loc(var, val)
1143       character(*) var, val
1144
1145 #ifdef WINIFL
1146       character(2000) line
1147       external ilen
1148
1149       open (196,file='env',status='old',readonly,shared)
1150       iread=0
1151 c      write(*,*)'looking for ',var
1152 10    read(196,*,err=11,end=11)line
1153       iread=index(line,var)
1154 c      write(*,*)iread,' ',var,' ',line
1155       if (iread.eq.0) go to 10
1156 c      write(*,*)'---> ',line
1157 11    continue
1158       if(iread.eq.0) then
1159 c       write(*,*)'CHUJ'
1160        val=''
1161       else
1162        iread=iread+ilen(var)+1
1163        read (line(iread:),*,err=12,end=12) val
1164 c       write(*,*)'OK: ',var,' = ',val
1165       endif
1166       close(196)
1167       return
1168 12    val=''
1169       close(196)
1170 #elif (defined CRAY)
1171       integer lennam,lenval,ierror
1172 c
1173 c        getenv using a POSIX call, useful on the T3D
1174 c        Sept 1996, comment out error check on advice of H. Pritchard
1175 c
1176       lennam = len(var)
1177       if(lennam.le.0) stop '--error calling getenv--'
1178       call pxfgetenv(var,lennam,val,lenval,ierror)
1179 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1180 #else
1181       call getenv(var,val)
1182 #endif
1183
1184       return
1185       end