source/unres/src_MD/parmread.F

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