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