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