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