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 c      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 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
 917      &  (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip(i),i=1,ntyp),
 918      &  (alp(i),i=1,ntyp)
 919 C For the GB potential convert sigma'**2 into chi'
 920       if (ipot.eq.4) then
 921         do i=1,ntyp
 922           chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
 923         enddo
 924       endif
 925       if (lprint) then
 926         write (iout,'(/a/)') 'Parameters of the BP potential:'
 927         write (iout,'(a/)') 'The epsilon array:'
 928         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 929         write (iout,'(/a)') 'One-body parameters:'
 930         write (iout,'(a,4x,4a)') 'residue','   sigma  ','s||/s_|_^2',
 931      &       '    chip  ','    alph  '
 932         write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
 933      &                     chip(i),alp(i),i=1,ntyp)
 934       endif
 935       goto 50
 936 C--------------------- GBV potential -----------------------------------
 937    40 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
 938      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
 939      &  (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
 940       if (lprint) then
 941         write (iout,'(/a/)') 'Parameters of the GBV potential:'
 942         write (iout,'(a/)') 'The epsilon array:'
 943         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 944         write (iout,'(/a)') 'One-body parameters:'
 945         write (iout,'(a,4x,5a)') 'residue','   sigma  ','    r0    ',
 946      &      's||/s_|_^2','    chip  ','    alph  '
 947         write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
 948      &           sigii(i),chip(i),alp(i),i=1,ntyp)
 949       endif
 950    50 continue
 951       close (isidep)
 952 C-----------------------------------------------------------------------
 953 C Calculate the "working" parameters of SC interactions.
 954       do i=2,ntyp
 955         do j=1,i-1
 956           eps(i,j)=eps(j,i)
 957         enddo
 958       enddo
 959       do i=1,ntyp
 960         do j=i,ntyp
 961           sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
 962           sigma(j,i)=sigma(i,j)
 963           rs0(i,j)=dwa16*sigma(i,j)
 964           rs0(j,i)=rs0(i,j)
 965         enddo
 966       enddo
 967       if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
 968      & 'Working parameters of the SC interactions:',
 969      & '     a    ','     b    ','   augm   ','  sigma ','   r0   ',
 970      & '  chi1   ','   chi2   '
 971       do i=1,ntyp
 972         do j=i,ntyp
 973           epsij=eps(i,j)
 974           if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
 975             rrij=sigma(i,j)
 976           else
 977             rrij=rr0(i)+rr0(j)
 978           endif
 979           r0(i,j)=rrij
 980           r0(j,i)=rrij
 981           rrij=rrij**expon
 982           epsij=eps(i,j)
 983           sigeps=dsign(1.0D0,epsij)
 984           epsij=dabs(epsij)
 985           aa(i,j)=epsij*rrij*rrij
 986           bb(i,j)=-sigeps*epsij*rrij
 987           aa(j,i)=aa(i,j)
 988           bb(j,i)=bb(i,j)
 989           if (ipot.gt.2) then
 990             sigt1sq=sigma0(i)**2
 991             sigt2sq=sigma0(j)**2
 992             sigii1=sigii(i)
 993             sigii2=sigii(j)
 994             ratsig1=sigt2sq/sigt1sq
 995             ratsig2=1.0D0/ratsig1
 996             chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
 997             if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
 998             rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
 999           else
1000             rsum_max=sigma(i,j)
1001           endif
1002 c         if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1003             sigmaii(i,j)=rsum_max
1004             sigmaii(j,i)=rsum_max
1005 c         else
1006 c           sigmaii(i,j)=r0(i,j)
1007 c           sigmaii(j,i)=r0(i,j)
1008 c         endif
1009 cd        write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1010           if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1011             r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1012             augm(i,j)=epsij*r_augm**(2*expon)
1013 c           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1014             augm(j,i)=augm(i,j)
1015           else
1016             augm(i,j)=0.0D0
1017             augm(j,i)=0.0D0
1018           endif
1019           if (lprint) then
1020             write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1021      &      restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1022      &      sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1023           endif
1024         enddo
1025       enddo
1026 #ifdef OLDSCP
1027 C
1028 C Define the SC-p interaction constants (hard-coded; old style)
1029 C
1030       do i=1,ntyp
1031 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1032 C helix formation)
1033 c       aad(i,1)=0.3D0*4.0D0**12
1034 C Following line for constants currently implemented
1035 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1036         aad(i,1)=1.5D0*4.0D0**12
1037 c       aad(i,1)=0.17D0*5.6D0**12
1038         aad(i,2)=aad(i,1)
1039 C "Soft" SC-p repulsion
1040         bad(i,1)=0.0D0
1041 C Following line for constants currently implemented
1042 c       aad(i,1)=0.3D0*4.0D0**6
1043 C "Hard" SC-p repulsion
1044         bad(i,1)=3.0D0*4.0D0**6
1045 c       bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1046         bad(i,2)=bad(i,1)
1047 c       aad(i,1)=0.0D0
1048 c       aad(i,2)=0.0D0
1049 c       bad(i,1)=1228.8D0
1050 c       bad(i,2)=1228.8D0
1051       enddo
1052 #else
1053 C
1054 C 8/9/01 Read the SC-p interaction constants from file
1055 C
1056       do i=1,ntyp
1057         read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1058       enddo
1059       do i=1,ntyp
1060         aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1061         aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1062         bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1063         bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1064       enddo
1065
1066       if (lprint) then
1067         write (iout,*) "Parameters of SC-p interactions:"
1068         do i=1,ntyp
1069           write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1070      &     eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1071         enddo
1072       endif
1073 #endif
1074 C
1075 C Define the constants of the disulfide bridge
1076 C
1077       ebr=-5.50D0
1078 c
1079 c Old arbitrary potential - commented out.
1080 c
1081 c      dbr= 4.20D0
1082 c      fbr= 3.30D0
1083 c
1084 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1085 c energy surface of diethyl disulfide.
1086 c A. Liwo and U. Kozlowska, 11/24/03
1087 c
1088       D0CM = 3.78d0
1089       AKCM = 15.1d0
1090       AKTH = 11.0d0
1091       AKCT = 12.0d0
1092       V1SS =-1.08d0
1093       V2SS = 7.61d0
1094       V3SS = 13.7d0
1095 c      akcm=0.0d0
1096 c      akth=0.0d0
1097 c      akct=0.0d0
1098 c      v1ss=0.0d0
1099 c      v2ss=0.0d0
1100 c      v3ss=0.0d0
1101
1102       if(me.eq.king) then
1103       write (iout,'(/a)') "Disulfide bridge parameters:"
1104       write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1105       write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1106       write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1107       write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1108      &  ' v3ss:',v3ss
1109       endif
1110       return
1111   111 write (iout,*) "Error reading bending energy parameters."
1112       goto 999
1113   112 write (iout,*) "Error reading rotamer energy parameters."
1114       goto 999
1115   113 write (iout,*) "Error reading torsional energy parameters."
1116       goto 999
1117   114 write (iout,*) "Error reading double torsional energy parameters."
1118       goto 999
1119   115 write (iout,*)
1120      &  "Error reading cumulant (multibody energy) parameters."
1121       goto 999
1122   116 write (iout,*) "Error reading electrostatic energy parameters."
1123       goto 999
1124   117 write (iout,*) "Error reading side chain interaction parameters."
1125       goto 999
1126   118 write (iout,*) "Error reading SCp interaction parameters."
1127       goto 999
1128   119 write (iout,*) "Error reading SCCOR parameters"
1129   999 continue
1130 #ifdef MPI
1131       call MPI_Finalize(Ierror)
1132 #endif
1133       stop
1134       return
1135       end
1136
1137
1138       subroutine getenv_loc(var, val)
1139       character(*) var, val
1140
1141 #ifdef WINIFL
1142       character(2000) line
1143       external ilen
1144
1145       open (196,file='env',status='old',readonly,shared)
1146       iread=0
1147 c      write(*,*)'looking for ',var
1148 10    read(196,*,err=11,end=11)line
1149       iread=index(line,var)
1150 c      write(*,*)iread,' ',var,' ',line
1151       if (iread.eq.0) go to 10
1152 c      write(*,*)'---> ',line
1153 11    continue
1154       if(iread.eq.0) then
1155 c       write(*,*)'CHUJ'
1156        val=''
1157       else
1158        iread=iread+ilen(var)+1
1159        read (line(iread:),*,err=12,end=12) val
1160 c       write(*,*)'OK: ',var,' = ',val
1161       endif
1162       close(196)
1163       return
1164 12    val=''
1165       close(196)
1166 #elif (defined CRAY)
1167       integer lennam,lenval,ierror
1168 c
1169 c        getenv using a POSIX call, useful on the T3D
1170 c        Sept 1996, comment out error check on advice of H. Pritchard
1171 c
1172       lennam = len(var)
1173       if(lennam.le.0) stop '--error calling getenv--'
1174       call pxfgetenv(var,lennam,val,lenval,ierror)
1175 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1176 #else
1177       call getenv(var,val)
1178 #endif
1179
1180       return
1181       end