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