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