source/unres/src_CSA_DiL/parmread.F

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