source/unres/src_MD/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       enddo
 324       write (2,*) "End reading THETA_PDB"
 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        itortyp(i)=-itortyp(-i)
 508       enddo
 509 c      write (iout,*) 'ntortyp',ntortyp
 510       do i=0,ntortyp-1
 511         do j=-ntortyp+1,ntortyp-1
 512           read (itorp,*,end=113,err=113) nterm(i,j,iblock),
 513      &          nlor(i,j,iblock)
 514           nterm(-i,-j,iblock)=nterm(i,j,iblock)
 515           nlor(-i,-j,iblock)=nlor(i,j,iblock)
 516           v0ij=0.0d0
 517           si=-1.0d0
 518           do k=1,nterm(i,j,iblock)
 519             read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),
 520      &      v2(k,i,j,iblock)
 521             v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
 522             v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
 523             v0ij=v0ij+si*v1(k,i,j,iblock)
 524             si=-si
 525 c         write(iout,*) i,j,k,iblock,nterm(i,j,iblock)
 526 c         write(iout,*) v1(k,-i,-j,iblock),v1(k,i,j,iblock),
 527 c      &v2(k,-i,-j,iblock),v2(k,i,j,iblock)
 528           enddo
 529           do k=1,nlor(i,j,iblock)
 530             read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),
 531      &        vlor2(k,i,j),vlor3(k,i,j)
 532             v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
 533           enddo
 534           v0(i,j,iblock)=v0ij
 535           v0(-i,-j,iblock)=v0ij
 536         enddo
 537       enddo
 538       enddo
 539       close (itorp)
 540       if (lprint) then
 541         write (iout,'(/a/)') 'Torsional constants:'
 542         do i=1,ntortyp
 543           do j=1,ntortyp
 544             write (iout,*) 'ityp',i,' jtyp',j
 545             write (iout,*) 'Fourier constants'
 546             do k=1,nterm(i,j,iblock)
 547               write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
 548      &        v2(k,i,j,iblock)
 549             enddo
 550             write (iout,*) 'Lorenz constants'
 551             do k=1,nlor(i,j,iblock)
 552               write (iout,'(3(1pe15.5))')
 553      &         vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
 554             enddo
 555           enddo
 556         enddo
 557       endif
 558 C
 559 C 6/23/01 Read parameters for double torsionals
 560 C
 561       do iblock=1,2
 562       do i=0,ntortyp-1
 563         do j=-ntortyp+1,ntortyp-1
 564           do k=-ntortyp+1,ntortyp-1
 565             read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
 566 c              write (iout,*) "OK onelett",
 567 c     &         i,j,k,t1,t2,t3
 568
 569             if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
 570      &        .or. t3.ne.toronelet(k)) then
 571               write (iout,*) "Error in double torsional parameter file",
 572      &         i,j,k,t1,t2,t3
 573 #ifdef MPI
 574               call MPI_Finalize(Ierror)
 575 #endif
 576                stop "Error in double torsional parameter file"
 577             endif
 578             read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),
 579      &         ntermd_2(i,j,k,iblock)
 580             ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
 581             ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
 582             read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,
 583      &         ntermd_1(i,j,k,iblock))
 584             read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,
 585      &         ntermd_1(i,j,k,iblock))
 586             read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,
 587      &         ntermd_1(i,j,k,iblock))
 588             read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,
 589      &         ntermd_1(i,j,k,iblock))
 590 C Martix of D parameters for one dimesional foureir series
 591             do l=1,ntermd_1(i,j,k,iblock)
 592              v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
 593              v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
 594              v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
 595              v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
 596 c            write(iout,*) "whcodze" ,
 597 c     & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
 598             enddo
 599             read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),
 600      &         v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
 601      &         v2s(m,l,i,j,k,iblock),
 602      &         m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
 603 C Martix of D parameters for two dimesional fourier series
 604             do l=1,ntermd_2(i,j,k,iblock)
 605              do m=1,l-1
 606              v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
 607              v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
 608              v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
 609              v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
 610              enddo!m
 611             enddo!l
 612           enddo!k
 613         enddo!j
 614       enddo!i
 615       enddo!iblock
 616       if (lprint) then
 617       write (iout,*)
 618       write (iout,*) 'Constants for double torsionals'
 619       do iblock=1,2
 620       do i=1,ntortyp
 621         do j=-ntortyp,ntortyp
 622           do k=-ntortyp,ntortyp
 623             write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
 624      &        ' nsingle',ntermd_1(i,j,k,iblock),
 625      &        ' ndouble',ntermd_2(i,j,k,iblock)
 626             write (iout,*)
 627             write (iout,*) 'Single angles:'
 628             do l=1,ntermd_1(i,j,k,iblock)
 629               write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
 630      &           v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
 631      &           v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),
 632      &           v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
 633             enddo
 634             write (iout,*)
 635             write (iout,*) 'Pairs of angles:'
 636             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
 637             do l=1,ntermd_2(i,j,k,iblock)
 638               write (iout,'(i5,20f10.5)')
 639      &         l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
 640             enddo
 641             write (iout,*)
 642             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
 643             do l=1,ntermd_2(i,j,k,iblock)
 644               write (iout,'(i5,20f10.5)')
 645      &         l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),
 646      &         (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
 647             enddo
 648             write (iout,*)
 649           enddo
 650         enddo
 651       enddo
 652       enddo
 653       endif
 654 #endif
 655 C Read of Side-chain backbone correlation parameters
 656 C Modified 11 May 2012 by Adasko
 657 CCC
 658 C
 659       read (isccor,*,end=113,err=113) nsccortyp
 660 #ifdef SCCORPDB
 661       read (isccor,*,end=113,err=113) (isccortyp(i),i=1,ntyp)
 662       do i=-ntyp,-1
 663         isccortyp(i)=-isccortyp(-i)
 664       enddo
 665       iscprol=isccortyp(20)
 666 c      write (iout,*) 'ntortyp',ntortyp
 667       maxinter=3
 668 cc maxinter is maximum interaction sites
 669       do l=1,maxinter
 670       do i=1,nsccortyp
 671         do j=1,nsccortyp
 672           read (isccor,*,end=113,err=113) nterm_sccor(i,j),nlor_sccor(i,j)
 673           v0ijsccor=0.0d0
 674           si=-1.0d0
 675
 676           do k=1,nterm_sccor(i,j)
 677             read (isccor,*,end=113,err=113) kk,v1sccor(k,l,i,j)
 678      &    ,v2sccor(k,l,i,j)
 679             if (j.eq.iscprol) then
 680              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
 681      &                        +v2sccor(k,l,i,j)*dsqrt(0.75d0)
 682              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
 683      &                        +v1sccor(k,l,i,j)*dsqrt(0.75d0)
 684              v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
 685              v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
 686              v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
 687              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
 688             else
 689              v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
 690              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
 691              v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
 692              v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
 693              v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
 694              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
 695             endif
 696             v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
 697             si=-si
 698           enddo
 699           do k=1,nlor_sccor(i,j)
 700             read (isccor,*,end=113,err=113) kk,vlor1sccor(k,i,j),
 701      &        vlor2sccor(k,i,j),vlor3sccor(k,i,j)
 702             v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
 703      &(1+vlor3sccor(k,i,j)**2)
 704           enddo
 705           v0sccor(i,j)=v0ijsccor
 706         enddo
 707       enddo
 708       enddo
 709       close (isccor)
 710 #else
 711       read (isccor,*,end=113,err=113) (isccortyp(i),i=1,ntyp)
 712 c      write (iout,*) 'ntortyp',ntortyp
 713       maxinter=3
 714 cc maxinter is maximum interaction sites
 715       do l=1,maxinter
 716       do i=1,nsccortyp
 717         do j=1,nsccortyp
 718           read (isccor,*,end=113,err=113)
 719       nterm_sccor(i,j),nlor_sccor(i,j)
 720           v0ijsccor=0.0d0
 721           si=-1.0d0
 722
 723           do k=1,nterm_sccor(i,j)
 724             read (isccor,*,end=113,err=113) kk,v1sccor(k,l,i,j)
 725      &    ,v2sccor(k,l,i,j)
 726             v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
 727             si=-si
 728           enddo
 729           do k=1,nlor_sccor(i,j)
 730             read (isccor,*,end=113,err=113) kk,vlor1sccor(k,i,j),
 731      &        vlor2sccor(k,i,j),vlor3sccor(k,i,j)
 732             v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
 733      &(1+vlor3sccor(k,i,j)**2)
 734           enddo
 735           v0sccor(i,j)=v0ijsccor
 736         enddo
 737       enddo
 738       enddo
 739       close (isccor)
 740
 741 #endif
 742       if (lprint) then
 743         write (iout,'(/a/)') 'Torsional constants:'
 744         do i=1,nsccortyp
 745           do j=1,nsccortyp
 746             write (iout,*) 'ityp',i,' jtyp',j
 747             write (iout,*) 'Fourier constants'
 748             do k=1,nterm_sccor(i,j)
 749       write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
 750             enddo
 751             write (iout,*) 'Lorenz constants'
 752             do k=1,nlor_sccor(i,j)
 753               write (iout,'(3(1pe15.5))')
 754      &         vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
 755             enddo
 756           enddo
 757         enddo
 758       endif
 759 C
 760 C
 761 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
 762 C         interaction energy of the Gly, Ala, and Pro prototypes.
 763 C
 764       if (lprint) then
 765         write (iout,*)
 766         write (iout,*) "Coefficients of the cumulants"
 767       endif
 768       read (ifourier,*) nloctyp
 769       do i=0,nloctyp-1
 770         read (ifourier,*,end=115,err=115)
 771         read (ifourier,*,end=115,err=115) (b(ii),ii=1,13)
 772         if (lprint) then
 773         write (iout,*) 'Type',i
 774         write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
 775         endif
 776         B1(1,i)  = b(3)
 777         B1(2,i)  = b(5)
 778         B1(1,-i) = b(3)
 779         B1(2,-i) = -b(5)
 780 c        b1(1,i)=0.0d0
 781 c        b1(2,i)=0.0d0
 782         B1tilde(1,i) = b(3)
 783         B1tilde(2,i) =-b(5)
 784         B1tilde(1,-i) =-b(3)
 785         B1tilde(2,-i) =b(5)
 786 c        b1tilde(1,i)=0.0d0
 787 c        b1tilde(2,i)=0.0d0
 788         B2(1,i)  = b(2)
 789         B2(2,i)  = b(4)
 790         B2(1,-i)  =b(2)
 791         B2(2,-i)  =-b(4)
 792
 793 c        b2(1,i)=0.0d0
 794 c        b2(2,i)=0.0d0
 795         CC(1,1,i)= b(7)
 796         CC(2,2,i)=-b(7)
 797         CC(2,1,i)= b(9)
 798         CC(1,2,i)= b(9)
 799         CC(1,1,-i)= b(7)
 800         CC(2,2,-i)=-b(7)
 801         CC(2,1,-i)=-b(9)
 802         CC(1,2,-i)=-b(9)
 803 c        CC(1,1,i)=0.0d0
 804 c        CC(2,2,i)=0.0d0
 805 c        CC(2,1,i)=0.0d0
 806 c        CC(1,2,i)=0.0d0
 807         Ctilde(1,1,i)=b(7)
 808         Ctilde(1,2,i)=b(9)
 809         Ctilde(2,1,i)=-b(9)
 810         Ctilde(2,2,i)=b(7)
 811         Ctilde(1,1,-i)=b(7)
 812         Ctilde(1,2,-i)=-b(9)
 813         Ctilde(2,1,-i)=b(9)
 814         Ctilde(2,2,-i)=b(7)
 815
 816 c        Ctilde(1,1,i)=0.0d0
 817 c        Ctilde(1,2,i)=0.0d0
 818 c        Ctilde(2,1,i)=0.0d0
 819 c        Ctilde(2,2,i)=0.0d0
 820         DD(1,1,i)= b(6)
 821         DD(2,2,i)=-b(6)
 822         DD(2,1,i)= b(8)
 823         DD(1,2,i)= b(8)
 824         DD(1,1,-i)= b(6)
 825         DD(2,2,-i)=-b(6)
 826         DD(2,1,-i)=-b(8)
 827         DD(1,2,-i)=-b(8)
 828 c        DD(1,1,i)=0.0d0
 829 c        DD(2,2,i)=0.0d0
 830 c        DD(2,1,i)=0.0d0
 831 c        DD(1,2,i)=0.0d0
 832         Dtilde(1,1,i)=b(6)
 833         Dtilde(1,2,i)=b(8)
 834         Dtilde(2,1,i)=-b(8)
 835         Dtilde(2,2,i)=b(6)
 836         Dtilde(1,1,-i)=b(6)
 837         Dtilde(1,2,-i)=-b(8)
 838         Dtilde(2,1,-i)=b(8)
 839         Dtilde(2,2,-i)=b(6)
 840
 841 c        Dtilde(1,1,i)=0.0d0
 842 c        Dtilde(1,2,i)=0.0d0
 843 c        Dtilde(2,1,i)=0.0d0
 844 c        Dtilde(2,2,i)=0.0d0
 845         EE(1,1,i)= b(10)+b(11)
 846         EE(2,2,i)=-b(10)+b(11)
 847         EE(2,1,i)= b(12)-b(13)
 848         EE(1,2,i)= b(12)+b(13)
 849         EE(1,1,-i)= b(10)+b(11)
 850         EE(2,2,-i)=-b(10)+b(11)
 851         EE(2,1,-i)=-b(12)+b(13)
 852         EE(1,2,-i)=-b(12)-b(13)
 853
 854 c        ee(1,1,i)=1.0d0
 855 c        ee(2,2,i)=1.0d0
 856 c        ee(2,1,i)=0.0d0
 857 c        ee(1,2,i)=0.0d0
 858 c        ee(2,1,i)=ee(1,2,i)
 859       enddo
 860       if (lprint) then
 861       do i=1,nloctyp
 862         write (iout,*) 'Type',i
 863         write (iout,*) 'B1'
 864         write(iout,*) B1(1,i),B1(2,i)
 865         write (iout,*) 'B2'
 866         write(iout,*) B2(1,i),B2(2,i)
 867         write (iout,*) 'CC'
 868         do j=1,2
 869           write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
 870         enddo
 871         write(iout,*) 'DD'
 872         do j=1,2
 873           write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
 874         enddo
 875         write(iout,*) 'EE'
 876         do j=1,2
 877           write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
 878         enddo
 879       enddo
 880       endif
 881 C
 882 C Read electrostatic-interaction parameters
 883 C
 884       if (lprint) then
 885         write (iout,*)
 886         write (iout,'(/a)') 'Electrostatic interaction constants:'
 887         write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
 888      &            'IT','JT','APP','BPP','AEL6','AEL3'
 889       endif
 890       read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
 891       read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
 892       read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
 893       read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
 894       close (ielep)
 895       do i=1,2
 896         do j=1,2
 897         rri=rpp(i,j)**6
 898         app (i,j)=epp(i,j)*rri*rri
 899         bpp (i,j)=-2.0D0*epp(i,j)*rri
 900         ael6(i,j)=elpp6(i,j)*4.2D0**6
 901         ael3(i,j)=elpp3(i,j)*4.2D0**3
 902         if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
 903      &                    ael6(i,j),ael3(i,j)
 904         enddo
 905       enddo
 906 C
 907 C Read side-chain interaction parameters.
 908 C
 909       read (isidep,*,end=117,err=117) ipot,expon
 910       if (ipot.lt.1 .or. ipot.gt.5) then
 911         write (iout,'(2a)') 'Error while reading SC interaction',
 912      &               'potential file - unknown potential type.'
 913 #ifdef MPI
 914         call MPI_Finalize(Ierror)
 915 #endif
 916         stop
 917       endif
 918       expon2=expon/2
 919       if(me.eq.king)
 920      & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
 921      & ', exponents are ',expon,2*expon
 922       goto (10,20,30,30,40) ipot
 923 C----------------------- LJ potential ---------------------------------
 924    10 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
 925      &   (sigma0(i),i=1,ntyp)
 926       if (lprint) then
 927         write (iout,'(/a/)') 'Parameters of the LJ potential:'
 928         write (iout,'(a/)') 'The epsilon array:'
 929         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 930         write (iout,'(/a)') 'One-body parameters:'
 931         write (iout,'(a,4x,a)') 'residue','sigma'
 932         write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
 933       endif
 934       goto 50
 935 C----------------------- LJK potential --------------------------------
 936    20 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
 937      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
 938       if (lprint) then
 939         write (iout,'(/a/)') 'Parameters of the LJK potential:'
 940         write (iout,'(a/)') 'The epsilon array:'
 941         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 942         write (iout,'(/a)') 'One-body parameters:'
 943         write (iout,'(a,4x,2a)') 'residue','   sigma  ','    r0    '
 944         write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
 945      &        i=1,ntyp)
 946       endif
 947       goto 50
 948 C---------------------- GB or BP potential -----------------------------
 949    30 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
 950      &  (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip(i),i=1,ntyp),
 951      &  (alp(i),i=1,ntyp)
 952 C For the GB potential convert sigma'**2 into chi'
 953       if (ipot.eq.4) then
 954         do i=1,ntyp
 955           chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
 956         enddo
 957       endif
 958       if (lprint) then
 959         write (iout,'(/a/)') 'Parameters of the BP potential:'
 960         write (iout,'(a/)') 'The epsilon array:'
 961         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 962         write (iout,'(/a)') 'One-body parameters:'
 963         write (iout,'(a,4x,4a)') 'residue','   sigma  ','s||/s_|_^2',
 964      &       '    chip  ','    alph  '
 965         write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
 966      &                     chip(i),alp(i),i=1,ntyp)
 967       endif
 968       goto 50
 969 C--------------------- GBV potential -----------------------------------
 970    40 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
 971      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
 972      &  (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
 973       if (lprint) then
 974         write (iout,'(/a/)') 'Parameters of the GBV potential:'
 975         write (iout,'(a/)') 'The epsilon array:'
 976         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 977         write (iout,'(/a)') 'One-body parameters:'
 978         write (iout,'(a,4x,5a)') 'residue','   sigma  ','    r0    ',
 979      &      's||/s_|_^2','    chip  ','    alph  '
 980         write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
 981      &           sigii(i),chip(i),alp(i),i=1,ntyp)
 982       endif
 983    50 continue
 984       close (isidep)
 985 C-----------------------------------------------------------------------
 986 C Calculate the "working" parameters of SC interactions.
 987       do i=2,ntyp
 988         do j=1,i-1
 989           eps(i,j)=eps(j,i)
 990         enddo
 991       enddo
 992       do i=1,ntyp
 993         do j=i,ntyp
 994           sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
 995           sigma(j,i)=sigma(i,j)
 996           rs0(i,j)=dwa16*sigma(i,j)
 997           rs0(j,i)=rs0(i,j)
 998         enddo
 999       enddo
1000       if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1001      & 'Working parameters of the SC interactions:',
1002      & '     a    ','     b    ','   augm   ','  sigma ','   r0   ',
1003      & '  chi1   ','   chi2   '
1004       do i=1,ntyp
1005         do j=i,ntyp
1006           epsij=eps(i,j)
1007           if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1008             rrij=sigma(i,j)
1009           else
1010             rrij=rr0(i)+rr0(j)
1011           endif
1012           r0(i,j)=rrij
1013           r0(j,i)=rrij
1014           rrij=rrij**expon
1015           epsij=eps(i,j)
1016           sigeps=dsign(1.0D0,epsij)
1017           epsij=dabs(epsij)
1018           aa(i,j)=epsij*rrij*rrij
1019           bb(i,j)=-sigeps*epsij*rrij
1020           aa(j,i)=aa(i,j)
1021           bb(j,i)=bb(i,j)
1022           if (ipot.gt.2) then
1023             sigt1sq=sigma0(i)**2
1024             sigt2sq=sigma0(j)**2
1025             sigii1=sigii(i)
1026             sigii2=sigii(j)
1027             ratsig1=sigt2sq/sigt1sq
1028             ratsig2=1.0D0/ratsig1
1029             chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1030             if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1031             rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1032           else
1033             rsum_max=sigma(i,j)
1034           endif
1035 c         if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1036             sigmaii(i,j)=rsum_max
1037             sigmaii(j,i)=rsum_max
1038 c         else
1039 c           sigmaii(i,j)=r0(i,j)
1040 c           sigmaii(j,i)=r0(i,j)
1041 c         endif
1042 cd        write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1043           if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1044             r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1045             augm(i,j)=epsij*r_augm**(2*expon)
1046 c           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1047             augm(j,i)=augm(i,j)
1048           else
1049             augm(i,j)=0.0D0
1050             augm(j,i)=0.0D0
1051           endif
1052           if (lprint) then
1053             write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1054      &      restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1055      &      sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1056           endif
1057         enddo
1058       enddo
1059 #ifdef OLDSCP
1060 C
1061 C Define the SC-p interaction constants (hard-coded; old style)
1062 C
1063       do i=1,20
1064 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1065 C helix formation)
1066 c       aad(i,1)=0.3D0*4.0D0**12
1067 C Following line for constants currently implemented
1068 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1069         aad(i,1)=1.5D0*4.0D0**12
1070 c       aad(i,1)=0.17D0*5.6D0**12
1071         aad(i,2)=aad(i,1)
1072 C "Soft" SC-p repulsion
1073         bad(i,1)=0.0D0
1074 C Following line for constants currently implemented
1075 c       aad(i,1)=0.3D0*4.0D0**6
1076 C "Hard" SC-p repulsion
1077         bad(i,1)=3.0D0*4.0D0**6
1078 c       bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1079         bad(i,2)=bad(i,1)
1080 c       aad(i,1)=0.0D0
1081 c       aad(i,2)=0.0D0
1082 c       bad(i,1)=1228.8D0
1083 c       bad(i,2)=1228.8D0
1084       enddo
1085 #else
1086 C
1087 C 8/9/01 Read the SC-p interaction constants from file
1088 C
1089       do i=1,ntyp
1090         read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1091       enddo
1092       do i=1,ntyp
1093         aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1094         aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1095         bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1096         bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1097       enddo
1098
1099       if (lprint) then
1100         write (iout,*) "Parameters of SC-p interactions:"
1101         do i=1,20
1102           write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1103      &     eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1104         enddo
1105       endif
1106 #endif
1107 C
1108 C Define the constants of the disulfide bridge
1109 C
1110       ebr=-5.50D0
1111 c
1112 c Old arbitrary potential - commented out.
1113 c
1114 c      dbr= 4.20D0
1115 c      fbr= 3.30D0
1116 c
1117 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1118 c energy surface of diethyl disulfide.
1119 c A. Liwo and U. Kozlowska, 11/24/03
1120 c
1121       D0CM = 3.78d0
1122       AKCM = 15.1d0
1123       AKTH = 11.0d0
1124       AKCT = 12.0d0
1125       V1SS =-1.08d0
1126       V2SS = 7.61d0
1127       V3SS = 13.7d0
1128 c      akcm=0.0d0
1129 c      akth=0.0d0
1130 c      akct=0.0d0
1131 c      v1ss=0.0d0
1132 c      v2ss=0.0d0
1133 c      v3ss=0.0d0
1134
1135       if(me.eq.king) then
1136       write (iout,'(/a)') "Disulfide bridge parameters:"
1137       write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1138       write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1139       write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1140       write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1141      &  ' v3ss:',v3ss
1142       endif
1143       return
1144   111 write (iout,*) "Error reading bending energy parameters."
1145       goto 999
1146   112 write (iout,*) "Error reading rotamer energy parameters."
1147       goto 999
1148   113 write (iout,*) "Error reading torsional energy parameters."
1149       goto 999
1150   114 write (iout,*) "Error reading double torsional energy parameters."
1151       goto 999
1152   115 write (iout,*)
1153      &  "Error reading cumulant (multibody energy) parameters."
1154       goto 999
1155   116 write (iout,*) "Error reading electrostatic energy parameters."
1156       goto 999
1157   117 write (iout,*) "Error reading side chain interaction parameters."
1158       goto 999
1159   118 write (iout,*) "Error reading SCp interaction parameters."
1160       goto 999
1161   119 write (iout,*) "Error reading SCCOR parameters"
1162   999 continue
1163 #ifdef MPI
1164       call MPI_Finalize(Ierror)
1165 #endif
1166       stop
1167       return
1168       end
1169
1170
1171       subroutine getenv_loc(var, val)
1172       character(*) var, val
1173
1174 #ifdef WINIFL
1175       character(2000) line
1176       external ilen
1177
1178       open (196,file='env',status='old',readonly,shared)
1179       iread=0
1180 c      write(*,*)'looking for ',var
1181 10    read(196,*,err=11,end=11)line
1182       iread=index(line,var)
1183 c      write(*,*)iread,' ',var,' ',line
1184       if (iread.eq.0) go to 10
1185 c      write(*,*)'---> ',line
1186 11    continue
1187       if(iread.eq.0) then
1188 c       write(*,*)'CHUJ'
1189        val=''
1190       else
1191        iread=iread+ilen(var)+1
1192        read (line(iread:),*,err=12,end=12) val
1193 c       write(*,*)'OK: ',var,' = ',val
1194       endif
1195       close(196)
1196       return
1197 12    val=''
1198       close(196)
1199 #elif (defined CRAY)
1200       integer lennam,lenval,ierror
1201 c
1202 c        getenv using a POSIX call, useful on the T3D
1203 c        Sept 1996, comment out error check on advice of H. Pritchard
1204 c
1205       lennam = len(var)
1206       if(lennam.le.0) stop '--error calling getenv--'
1207       call pxfgetenv(var,lennam,val,lenval,ierror)
1208 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1209 #else
1210       call getenv(var,val)
1211 #endif
1212
1213       return
1214       end