source/unres/src_MD-NEWSC/parmread_v3ok1.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
 557 C 5/21/07 (AL) Read coefficients of the backbone-local sidechain-local
 558 C         correlation energies.
 559 C
 560       read (isccor,*,end=119,err=119) nterm_sccor
 561       do i=1,20
 562         do j=1,20
 563           read (isccor,'(a)')
 564           do k=1,nterm_sccor
 565             read (isccor,*,end=119,err=119) kk,v1sccor(k,i,j),
 566      &        v2sccor(k,i,j)
 567           enddo
 568         enddo
 569       enddo
 570       close (isccor)
 571       if (lprint) then
 572         write (iout,'(/a/)') 'Torsional constants of SCCORR:'
 573         do i=1,20
 574           do j=1,20
 575             write (iout,*) 'ityp',i,' jtyp',j
 576             do k=1,nterm_sccor
 577               write (iout,'(2(1pe15.5))') v1sccor(k,i,j),v2sccor(k,i,j)
 578             enddo
 579           enddo
 580         enddo
 581       endif
 582 C
 583 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
 584 C         interaction energy of the Gly, Ala, and Pro prototypes.
 585 C
 586       if (lprint) then
 587         write (iout,*)
 588         write (iout,*) "Coefficients of the cumulants"
 589       endif
 590       read (ifourier,*) nloctyp
 591       do i=1,nloctyp
 592         read (ifourier,*,end=115,err=115)
 593         read (ifourier,*,end=115,err=115) (b(ii),ii=1,13)
 594         if (lprint) then
 595         write (iout,*) 'Type',i
 596         write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
 597         endif
 598         B1(1,i)  = b(3)
 599         B1(2,i)  = b(5)
 600 c        b1(1,i)=0.0d0
 601 c        b1(2,i)=0.0d0
 602         B1tilde(1,i) = b(3)
 603         B1tilde(2,i) =-b(5)
 604 c        b1tilde(1,i)=0.0d0
 605 c        b1tilde(2,i)=0.0d0
 606         B2(1,i)  = b(2)
 607         B2(2,i)  = b(4)
 608 c        b2(1,i)=0.0d0
 609 c        b2(2,i)=0.0d0
 610         CC(1,1,i)= b(7)
 611         CC(2,2,i)=-b(7)
 612         CC(2,1,i)= b(9)
 613         CC(1,2,i)= b(9)
 614 c        CC(1,1,i)=0.0d0
 615 c        CC(2,2,i)=0.0d0
 616 c        CC(2,1,i)=0.0d0
 617 c        CC(1,2,i)=0.0d0
 618         Ctilde(1,1,i)=b(7)
 619         Ctilde(1,2,i)=b(9)
 620         Ctilde(2,1,i)=-b(9)
 621         Ctilde(2,2,i)=b(7)
 622 c        Ctilde(1,1,i)=0.0d0
 623 c        Ctilde(1,2,i)=0.0d0
 624 c        Ctilde(2,1,i)=0.0d0
 625 c        Ctilde(2,2,i)=0.0d0
 626         DD(1,1,i)= b(6)
 627         DD(2,2,i)=-b(6)
 628         DD(2,1,i)= b(8)
 629         DD(1,2,i)= b(8)
 630 c        DD(1,1,i)=0.0d0
 631 c        DD(2,2,i)=0.0d0
 632 c        DD(2,1,i)=0.0d0
 633 c        DD(1,2,i)=0.0d0
 634         Dtilde(1,1,i)=b(6)
 635         Dtilde(1,2,i)=b(8)
 636         Dtilde(2,1,i)=-b(8)
 637         Dtilde(2,2,i)=b(6)
 638 c        Dtilde(1,1,i)=0.0d0
 639 c        Dtilde(1,2,i)=0.0d0
 640 c        Dtilde(2,1,i)=0.0d0
 641 c        Dtilde(2,2,i)=0.0d0
 642         EE(1,1,i)= b(10)+b(11)
 643         EE(2,2,i)=-b(10)+b(11)
 644         EE(2,1,i)= b(12)-b(13)
 645         EE(1,2,i)= b(12)+b(13)
 646 c        ee(1,1,i)=1.0d0
 647 c        ee(2,2,i)=1.0d0
 648 c        ee(2,1,i)=0.0d0
 649 c        ee(1,2,i)=0.0d0
 650 c        ee(2,1,i)=ee(1,2,i)
 651       enddo
 652       if (lprint) then
 653       do i=1,nloctyp
 654         write (iout,*) 'Type',i
 655         write (iout,*) 'B1'
 656         write(iout,*) B1(1,i),B1(2,i)
 657         write (iout,*) 'B2'
 658         write(iout,*) B2(1,i),B2(2,i)
 659         write (iout,*) 'CC'
 660         do j=1,2
 661           write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
 662         enddo
 663         write(iout,*) 'DD'
 664         do j=1,2
 665           write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
 666         enddo
 667         write(iout,*) 'EE'
 668         do j=1,2
 669           write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
 670         enddo
 671       enddo
 672       endif
 673 C
 674 C Read electrostatic-interaction parameters
 675 C
 676       if (lprint) then
 677         write (iout,*)
 678         write (iout,'(/a)') 'Electrostatic interaction constants:'
 679         write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
 680      &            'IT','JT','APP','BPP','AEL6','AEL3'
 681       endif
 682       read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
 683       read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
 684       read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
 685       read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
 686       close (ielep)
 687       do i=1,2
 688         do j=1,2
 689         rri=rpp(i,j)**6
 690         app (i,j)=epp(i,j)*rri*rri
 691         bpp (i,j)=-2.0D0*epp(i,j)*rri
 692         ael6(i,j)=elpp6(i,j)*4.2D0**6
 693         ael3(i,j)=elpp3(i,j)*4.2D0**3
 694         if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
 695      &                    ael6(i,j),ael3(i,j)
 696         enddo
 697       enddo
 698 C
 699 C Read side-chain interaction parameters.
 700 C
 701       read (isidep,*,end=117,err=117) ipot,expon
 702       if (ipot.lt.1 .or. ipot.gt.6) then
 703         write (iout,'(2a)') 'Error while reading SC interaction',
 704      &               'potential file - unknown potential type.'
 705 #ifdef MPI
 706         call MPI_Finalize(Ierror)
 707 #endif
 708         stop
 709       endif
 710       expon2=expon/2
 711       if(me.eq.king)
 712      & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
 713      & ', exponents are ',expon,2*expon
 714       goto (10,20,30,30,40,50) ipot
 715 C----------------------- LJ potential ---------------------------------
 716    10 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
 717      &   (sigma0(i),i=1,ntyp)
 718       if (lprint) then
 719         write (iout,'(/a/)') 'Parameters of the LJ potential:'
 720         write (iout,'(a/)') 'The epsilon array:'
 721         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 722         write (iout,'(/a)') 'One-body parameters:'
 723         write (iout,'(a,4x,a)') 'residue','sigma'
 724         write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
 725       endif
 726       goto 60
 727 C----------------------- LJK potential --------------------------------
 728    20 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
 729      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
 730       if (lprint) then
 731         write (iout,'(/a/)') 'Parameters of the LJK potential:'
 732         write (iout,'(a/)') 'The epsilon array:'
 733         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 734         write (iout,'(/a)') 'One-body parameters:'
 735         write (iout,'(a,4x,2a)') 'residue','   sigma  ','    r0    '
 736         write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
 737      &        i=1,ntyp)
 738       endif
 739       goto 60
 740 C---------------------- GB or BP potential -----------------------------
 741    30 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
 742      &  (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip(i),i=1,ntyp),
 743      &  (alp(i),i=1,ntyp)
 744 C For the GB potential convert sigma'**2 into chi'
 745       if (ipot.eq.4) then
 746         do i=1,ntyp
 747           chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
 748         enddo
 749       endif
 750       if (lprint) then
 751         write (iout,'(/a/)') 'Parameters of the BP potential:'
 752         write (iout,'(a/)') 'The epsilon array:'
 753         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 754         write (iout,'(/a)') 'One-body parameters:'
 755         write (iout,'(a,4x,4a)') 'residue','   sigma  ','s||/s_|_^2',
 756      &       '    chip  ','    alph  '
 757         write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
 758      &                     chip(i),alp(i),i=1,ntyp)
 759       endif
 760       goto 60
 761 C--------------------- GBV potential -----------------------------------
 762    40 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
 763      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
 764      &  (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
 765       if (lprint) then
 766         write (iout,'(/a/)') 'Parameters of the GBV potential:'
 767         write (iout,'(a/)') 'The epsilon array:'
 768         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 769         write (iout,'(/a)') 'One-body parameters:'
 770         write (iout,'(a,4x,5a)') 'residue','   sigma  ','    r0    ',
 771      &      's||/s_|_^2','    chip  ','    alph  '
 772         write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
 773      &           sigii(i),chip(i),alp(i),i=1,ntyp)
 774       endif
 775       goto 60
 776
 777 C--------------------- Momo potential -----------------------------------
 778
 779    50 continue
 780
 781       read (isidep,*) (icharge(i),i=1,ntyp)
 782 c      write (2,*) "icharge",(icharge(i),i=1,ntyp)
 783       do i=1,ntyp
 784        do j=1,i
 785 c!        write (*,*) "Im in ", i, " ", j
 786         read(isidep,*)
 787      &  eps(i,j),sigma(i,j),chi(i,j),chi(j,i),chipp(i,j),chipp(j,i),
 788      &  (alphasur(k,i,j),k=1,4),sigmap(i,j),sigmap(j,i),
 789      &  chis(i,j),chis(j,i),
 790      &  nstate(i,j),(wstate(k,i,j),k=1,4),
 791      &  ((dhead(l,k,i,j),k=1,2),l=1,2),dtail(i,j),dtail(j,i),
 792      &  epshead(i,j),sig0head(i,j),
 793      &  rborn(i,j),rborn(j,i),(wqdip(k,i,j),k=1,2),wquad(i,j),
 794      &  alphapol(i,j),alphapol(j,i),
 795      &  (alphiso(k,i,j),k=1,4),sigiso(i,j),sigiso(j,i),epsintab(i,j)
 796        END DO
 797       END DO
 798
 799 c!      write (*,*) "Parameters read in"
 800
 801 c! THIS LOOP FILLS PARAMETERS FOR PAIRS OF AA's NOT EXPLICITLY
 802 c! DEFINED IN SCPARM.MOMO. IT DOES SO BY TAKING THEM FROM SYMMETRIC
 803 c! PAIR, FOR EG. IF ARG-HIS IS BLANK, IT WILL TRY TO TAKE PARAMETERS
 804 c! FROM HIS-ARG.
 805 c!
 806 c! DISABLE IT AT >>YOUR OWN PERIL<<
 807 c!
 808       DO i = 1, ntyp
 809        DO j = i+1, ntyp
 810         eps(i,j) = eps(j,i)
 811         sigma(i,j) = sigma(j,i)
 812         nstate(i,j) = nstate(j,i)
 813
 814         DO k = 1, 4
 815          alphasur(k,i,j) = alphasur(k,j,i)
 816          wstate(k,i,j) = wstate(k,j,i)
 817          alphiso(k,i,j) = alphiso(k,j,i)
 818         END DO
 819
 820         DO k = 1, 2
 821          DO l = 1, 2
 822           dhead(k,l,i,j) = dhead(k,l,j,i)
 823          END DO
 824         END DO
 825
 826         epshead(i,j) = epshead(j,i)
 827         sig0head(i,j) = sig0head(j,i)
 828
 829         DO k = 1, 2
 830          wqdip(k,i,j) = wqdip(k,j,i)
 831         END DO
 832
 833         wquad(i,j) = wquad(j,i)
 834         epsintab(i,j) = epsintab(j,i)
 835
 836        END DO
 837       END DO
 838
 839 c!      IF (1.eq.0) THEN
 840 c!      write (*,*) ""
 841 c!      write (*,*) ""
 842 c!      write (*,'(f17.15)') ((eps(i,j),j=i,ntyp), i=1, ntyp)
 843 c!      write (*,'(f17.15)') (sigma0(i),i=1,ntyp)
 844 c!      write (*,'(f17.15)') (sigii(i),i=1,ntyp)
 845 c!      write (*,'(f17.15)') (chip(i),i=1,ntyp)
 846 c!      write (*,'(f17.15)') (alp(i),i=1,ntyp)
 847 c!      write (*,*) ""
 848 c!      write (*,*) ""
 849 c!      END IF
 850
 851 c! This is older parameter-filling loop
 852 c!      DO i = 1, ntyp
 853 c!       DO j = 1, ntyp
 854 c
 855 c        IF ((sig0head(i,j).EQ.0.0d0).AND.(sig0head(j,i).NE.0.0d0)) THEN
 856 c         sig0head(i,j) = sig0head(j,i)
 857 c        ELSE IF ((sig0head(j,i).EQ.0.0d0).AND.(sig0head(i,j).NE.0.0d0))
 858 c     &  THEN
 859 c         sig0head(j,i) = sig0head(i,j)
 860 c        END IF
 861 c
 862 c        IF ((epshead(i,j).EQ.0.0d0).AND.(epshead(j,i).NE.0.0d0)) THEN
 863 c         epshead(i,j) = epshead(j,i)
 864 c        ELSE IF ((epshead(j,i).EQ.0.0d0).AND.(epshead(i,j).NE.0.0d0))
 865 c     &  THEN
 866 c         epshead(j,i) = epshead(i,j)
 867 c        END IF
 868 c
 869 c        IF ((wquad(i,j).EQ.0.0d0).AND.(wquad(j,i).NE.0.0d0)) THEN
 870 c         wquad(i,j) = wquad(j,i)
 871 c        ELSE IF ((wquad(j,i).EQ.0.0d0).AND.(wquad(i,j).NE.0.0d0))
 872 c     &  THEN
 873 c         wquad(j,i) = wquad(i,j)
 874 c        END IF
 875 c
 876 c        IF ((sigiso(i,j).EQ.0.0d0).AND.sigiso(j,i).NE.0.0d0)) THEN
 877 c         sigiiso(i,j) = sigiso(j,i)
 878 c        ELSE IF ((epshead(j,i).EQ.0.0d0).AND.(epshead(i,j).NE.0.0d0))
 879 c     &  THEN
 880 c         sigiso(j,i) = sigiso(i,j)
 881 c        END IF
 882 c!        DO k = 1, 4
 883 c!         IF ((wstate(k,i,j).EQ.0.0d0).AND.
 884 c!     &   (wstate(k,j,i).NE.0.0d0)) THEN
 885 c!          wstate(k,i,j) = wstate(k,j,i)
 886 c!         ELSE IF ((wstate(k,j,i).EQ.0.0d0).AND.
 887 c!     &   (wstate(k,i,j).NE.0.0d0)) THEN
 888 c!          wstate(k,j,i) = wstate(k,i,j)
 889 c!         END IF
 890 c!        END DO
 891 c!       END DO
 892 c!      END DO
 893 c! THE LOOP MAKING MOMO_TABLES SYMMETRIC ENDS HERE
 894
 895       if (.not.lprint) goto 70
 896       write (iout,'(a)')
 897      & "Parameters of the new physics-based SC-SC interaction potential"
 898       write (iout,'(/7a)') 'Residues','     epsGB','       rGB',
 899      &  '    chi1GB','    chi2GB','   chip1GB','   chip2GB'
 900       do i=1,ntyp
 901         do j=1,i
 902           write (iout,'(2(a3,1x),1pe10.3,5(0pf10.3))')
 903      &      restyp(i),restyp(j),eps(i,j),sigma(i,j),chi(i,j),chi(j,i),
 904      &       chipp(i,j),chipp(j,i)
 905         enddo
 906       enddo
 907       write (iout,'(/9a)') 'Residues',' alphasur1',' alphasur2',
 908      &  ' alphasur3',' alphasur4','   sigmap1','   sigmap2',
 909      &  '     chis1','     chis2'
 910       do i=1,ntyp
 911         do j=1,i
 912           write (iout,'(2(a3,1x),8(0pf10.3))')
 913      &      restyp(i),restyp(j),(alphasur(k,i,j),k=1,4),
 914      &       sigmap(i,j),sigmap(j,i),chis(i,j),chis(j,i)
 915         enddo
 916       enddo
 917       write (iout,'(/14a)') 'Residues',' nst ','  wst1',
 918      &  '    wst2','    wst3','    wst4','   dh11','   dh21',
 919      &  '   dh12','   dh22','    dt1','    dt2','    epsh1',
 920      &  '   sigh'
 921       do i=1,ntyp
 922         do j=1,i
 923           write (iout,'(2(a3,1x),i3,4f8.3,6f7.2,f9.5,f7.2)')
 924      &      restyp(i),restyp(j),nstate(i,j),(wstate(k,i,j),k=1,4),
 925      &       ((dhead(l,k,i,j),l=1,2),k=1,2),dtail(i,j),dtail(j,i),
 926      &       epshead(i,j),sig0head(i,j)
 927         enddo
 928       enddo
 929       write (iout,'(/12a)') 'Residues',' ch1',' ch2',
 930      &  '    rborn1','    rborn2','    wqdip1','    wqdip2',
 931      &  '     wquad'
 932       do i=1,ntyp
 933         do j=1,i
 934           write (iout,'(2(a3,1x),2i4,5f10.3)')
 935      &      restyp(i),restyp(j),icharge(i),icharge(j),
 936      &      rborn(i,j),rborn(j,i),(wqdip(k,i,j),k=1,2),wquad(i,j)
 937         enddo
 938       enddo
 939       write (iout,'(/12a)') 'Residues',
 940      &  '  alphpol1',
 941      &  '  alphpol2','  alphiso1','   alpiso2',
 942      &  '   alpiso3','   alpiso4','   sigiso1','   sigiso2',
 943      &  '     epsin'
 944       do i=1,ntyp
 945         do j=1,i
 946           write (iout,'(2(a3,1x),11f10.3)')
 947      &      restyp(i),restyp(j),alphapol(i,j),alphapol(j,i),
 948      &      (alphiso(k,i,j),k=1,4),sigiso(i,j),sigiso(j,i),epsintab(i,j)
 949         enddo
 950       enddo
 951       goto 70
 952
 953 C For the GB potential convert sigma'**2 into chi'
 954       DO i=1,ntyp
 955        chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
 956       END DO
 957       IF (lprint) THEN
 958        write (iout,'(/a/)') 'Parameters of the BP potential:'
 959        write (iout,'(a/)') 'The epsilon array:'
 960        CALL printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 961        write (iout,'(/a)') 'One-body parameters:'
 962        write (iout,'(a,4x,4a)') 'residue','   sigma  ','s||/s_|_^2',
 963      &       '    chip  ','    alph  '
 964        write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
 965      &                     chip(i),alp(i),i=1,ntyp)
 966       END IF
 967       GOTO 60
 968    60 CONTINUE
 969       close (isidep)
 970 C-----------------------------------------------------------------------
 971 C Calculate the "working" parameters of SC interactions.
 972       do i=2,ntyp
 973         do j=1,i-1
 974           eps(i,j)=eps(j,i)
 975           do k=1,4
 976             alphasur(k,j,i)=alphasur(k,i,j)
 977             alphiso(k,j,i)=alphiso(k,i,j)
 978             wstate(k,j,i)=wstate(k,i,j)
 979           enddo
 980           do k=1,2
 981             wqdip(k,j,i)=wqdip(k,i,j)
 982           enddo
 983           do k=1,2
 984             do l=1,2
 985               dhead(l,k,j,i)=dhead(l,k,i,j)
 986             enddo
 987           enddo
 988         enddo
 989       enddo
 990       IF (ipot.LT.6) THEN
 991        do i=1,ntyp
 992         do j=i,ntyp
 993          sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
 994          sigma(j,i)=sigma(i,j)
 995          rs0(i,j)=dwa16*sigma(i,j)
 996          rs0(j,i)=rs0(i,j)
 997         enddo
 998        enddo
 999       END IF
1000
1001    70 continue
1002
1003       if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1004      & 'Working parameters of the SC interactions:',
1005      & '     a    ','     b    ','   augm   ','  sigma ','   r0   ',
1006      & '  chi1   ','   chi2   '
1007       do i=1,ntyp
1008         do j=i,ntyp
1009           epsij=eps(i,j)
1010           IF (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4 .or. ipot.eq.6)
1011      &    THEN
1012             rrij=sigma(i,j)
1013           ELSE
1014             rrij=rr0(i)+rr0(j)
1015           END IF
1016           r0(i,j)=rrij
1017           r0(j,i)=rrij
1018           rrij=rrij**expon
1019           epsij=eps(i,j)
1020           sigeps=dsign(1.0D0,epsij)
1021           epsij=dabs(epsij)
1022           aa(i,j)=epsij*rrij*rrij
1023           bb(i,j)=-sigeps*epsij*rrij
1024           aa(j,i)=aa(i,j)
1025           bb(j,i)=bb(i,j)
1026           IF ((ipot.GT.2).AND.(ipot.LT.6)) THEN
1027             sigt1sq  = sigma0(i)**2
1028             sigt2sq  = sigma0(j)**2
1029             sigii1   = sigii(i)
1030             sigii2   = sigii(j)
1031             ratsig1  = sigt2sq/sigt1sq
1032             ratsig2  = 1.0D0/ratsig1
1033             chi(i,j) = (sigii1-1.0D0) / (sigii1+ratsig1)
1034             IF (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1035             rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1036           ELSE
1037             rsum_max=sigma(i,j)
1038           END IF
1039 c         if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1040             sigmaii(i,j)=rsum_max
1041             sigmaii(j,i)=rsum_max
1042 c         else
1043 c           sigmaii(i,j)=r0(i,j)
1044 c           sigmaii(j,i)=r0(i,j)
1045 c         endif
1046 cd        write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1047           if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1048             r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1049             augm(i,j)=epsij*r_augm**(2*expon)
1050 c           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1051             augm(j,i)=augm(i,j)
1052           else
1053             augm(i,j)=0.0D0
1054             augm(j,i)=0.0D0
1055           endif
1056
1057           if (lprint) then
1058             if (ipot.lt.6) then
1059             write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1060      &      restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1061      &      sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1062             else
1063             write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3),2i3,10f8.4,
1064      &       i3,40f10.4)')
1065      &      restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1066      &      sigma(i,j),r0(i,j),chi(i,j),chi(j,i),
1067      &      icharge(i),icharge(j),chipp(i,j),chipp(j,i),
1068      &     (alphasur(k,i,j),k=1,4),sigmap(i,j),sigmap(j,i),
1069      &     chis(i,j),chis(j,i),
1070      &     nstate(i,j),(wstate(k,i,j),k=1,4),
1071      &     ((dhead(l,k,i,j),l=1,2),k=1,2),dtail(i,j),dtail(j,i),
1072      &     epshead(i,j),sig0head(i,j),
1073      &     rborn(i,j),(wqdip(k,i,j),k=1,2),wquad(i,j),
1074      &     alphapol(i,j),alphapol(j,i),
1075      &     (alphiso(k,i,j),k=1,4),sigiso(i,j)
1076
1077             endif
1078                 endif
1079         enddo
1080       enddo
1081
1082
1083
1084 #ifdef OLDSCP
1085 C
1086 C Define the SC-p interaction constants (hard-coded; old style)
1087 C
1088       do i=1,20
1089 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1090 C helix formation)
1091 c       aad(i,1)=0.3D0*4.0D0**12
1092 C Following line for constants currently implemented
1093 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1094         aad(i,1)=1.5D0*4.0D0**12
1095 c       aad(i,1)=0.17D0*5.6D0**12
1096         aad(i,2)=aad(i,1)
1097 C "Soft" SC-p repulsion
1098         bad(i,1)=0.0D0
1099 C Following line for constants currently implemented
1100 c       aad(i,1)=0.3D0*4.0D0**6
1101 C "Hard" SC-p repulsion
1102         bad(i,1)=3.0D0*4.0D0**6
1103 c       bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1104         bad(i,2)=bad(i,1)
1105 c       aad(i,1)=0.0D0
1106 c       aad(i,2)=0.0D0
1107 c       bad(i,1)=1228.8D0
1108 c       bad(i,2)=1228.8D0
1109       enddo
1110 #else
1111 C
1112 C 8/9/01 Read the SC-p interaction constants from file
1113 C
1114       do i=1,ntyp
1115         read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1116       enddo
1117       do i=1,ntyp
1118         aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1119         aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1120         bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1121         bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1122       enddo
1123
1124       if (lprint) then
1125         write (iout,'(/a)') "Parameters of SC-p interactions:"
1126         do i=1,20
1127           write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1128      &     eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1129         enddo
1130       endif
1131 #endif
1132 C
1133 C Define the constants of the disulfide bridge
1134 C
1135       ebr=-5.50D0
1136 c
1137 c Old arbitrary potential - commented out.
1138 c
1139 c      dbr= 4.20D0
1140 c      fbr= 3.30D0
1141 c
1142 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1143 c energy surface of diethyl disulfide.
1144 c A. Liwo and U. Kozlowska, 11/24/03
1145 c
1146       D0CM = 3.78d0
1147       AKCM = 15.1d0
1148       AKTH = 11.0d0
1149       AKCT = 12.0d0
1150       V1SS =-1.08d0
1151       V2SS = 7.61d0
1152       V3SS = 13.7d0
1153 c      akcm=0.0d0
1154 c      akth=0.0d0
1155 c      akct=0.0d0
1156 c      v1ss=0.0d0
1157 c      v2ss=0.0d0
1158 c      v3ss=0.0d0
1159
1160 c!      print *, " "
1161 c!      print *, "END OF PARMREAD"
1162 c!      print *, "eps = ", eps(9,9), "sigma = ", sigma(9,9)
1163 c!      print *, "chi1 = ", chi(9,9), "chi2 = ", chi(9,9)
1164 c!      print *, "chip1 = ", chipp(9,9), "chip2 = ", chipp(9,9)
1165 c!      print *, "sig1 = ", sigmap(9,9), "sig2 = ", sigmap(9,9)
1166 c!      print *, "chis1 = ", chis(9,9)," chis2 = ", chis(9,9)
1167 c!      print *, "END OF PARMREAD"
1168 c!      print *, " "
1169
1170
1171       if(me.eq.king) then
1172       write (iout,'(/a)') "Disulfide bridge parameters:"
1173       write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1174       write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1175       write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1176       write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1177      &  ' v3ss:',v3ss
1178       endif
1179       return
1180   111 write (iout,*) "Error reading bending energy parameters."
1181       goto 999
1182   112 write (iout,*) "Error reading rotamer energy parameters."
1183       goto 999
1184   113 write (iout,*) "Error reading torsional energy parameters."
1185       goto 999
1186   114 write (iout,*) "Error reading double torsional energy parameters."
1187       goto 999
1188   115 write (iout,*)
1189      &  "Error reading cumulant (multibody energy) parameters."
1190       goto 999
1191   116 write (iout,*) "Error reading electrostatic energy parameters."
1192       goto 999
1193   117 write (iout,*) "Error reading side chain interaction parameters."
1194       goto 999
1195   118 write (iout,*) "Error reading SCp interaction parameters."
1196       goto 999
1197   119 write (iout,*) "Error reading SCCOR parameters"
1198   999 continue
1199 #ifdef MPI
1200       call MPI_Finalize(Ierror)
1201 #endif
1202       stop
1203       return
1204       end
1205
1206
1207       subroutine getenv_loc(var, val)
1208       character(*) var, val
1209
1210 #ifdef WINIFL
1211       character(2000) line
1212       external ilen
1213
1214       open (196,file='env',status='old',readonly,shared)
1215       iread=0
1216 c      write(*,*)'looking for ',var
1217 10    read(196,*,err=11,end=11)line
1218       iread=index(line,var)
1219 c      write(*,*)iread,' ',var,' ',line
1220       if (iread.eq.0) go to 10
1221 c      write(*,*)'---> ',line
1222 11    continue
1223       if(iread.eq.0) then
1224 c       write(*,*)'CHUJ'
1225        val=''
1226       else
1227        iread=iread+ilen(var)+1
1228        read (line(iread:),*,err=12,end=12) val
1229 c       write(*,*)'OK: ',var,' = ',val
1230       endif
1231       close(196)
1232       return
1233 12    val=''
1234       close(196)
1235 #elif (defined CRAY)
1236       integer lennam,lenval,ierror
1237 c
1238 c        getenv using a POSIX call, useful on the T3D
1239 c        Sept 1996, comment out error check on advice of H. Pritchard
1240 c
1241       lennam = len(var)
1242       if(lennam.le.0) stop '--error calling getenv--'
1243       call pxfgetenv(var,lennam,val,lenval,ierror)
1244 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1245 #else
1246       call getenv(var,val)
1247 #endif
1248
1249       return
1250       end