source/unres/src_MD-M/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       print *,i
  77         read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),
  78      &   j=1,nbondterm(i)),msc(i),isc(i),restok(i)
  79         dsc(i) = vbldsc0(1,i)
  80         if (i.eq.10) then
  81           dsc_inv(i)=0.0D0
  82         else
  83           dsc_inv(i)=1.0D0/dsc(i)
  84         endif
  85       enddo
  86 #endif
  87       if (lprint) then
  88         write(iout,'(/a/)')"Dynamic constants of the interaction sites:"
  89         write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K','A0','mass',
  90      &   'inertia','Pstok'
  91         write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0,mp,ip,pstok
  92         do i=1,ntyp
  93           write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),
  94      &      vbldsc0(1,i),aksc(1,i),abond0(1,i),msc(i),isc(i),restok(i)
  95           do j=2,nbondterm(i)
  96             write (iout,'(13x,3f10.5)')
  97      &        vbldsc0(j,i),aksc(j,i),abond0(j,i)
  98           enddo
  99         enddo
 100       endif
 101 #ifdef CRYST_THETA
 102 C
 103 C Read the parameters of the probability distribution/energy expression
 104 C of the virtual-bond valence angles theta
 105 C
 106       do i=1,ntyp
 107         read (ithep,*,err=111,end=111) a0thet(i),(athet(j,i,1,1),j=1,2),
 108      &    (bthet(j,i,1,1),j=1,2)
 109         read (ithep,*,err=111,end=111) (polthet(j,i),j=0,3)
 110         read (ithep,*,err=111,end=111) (gthet(j,i),j=1,3)
 111         read (ithep,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
 112         sigc0(i)=sigc0(i)**2
 113       enddo
 114       do i=1,ntyp
 115       athet(1,i,1,-1)=athet(1,i,1,1)
 116       athet(2,i,1,-1)=athet(2,i,1,1)
 117       bthet(1,i,1,-1)=-bthet(1,i,1,1)
 118       bthet(2,i,1,-1)=-bthet(2,i,1,1)
 119       athet(1,i,-1,1)=-athet(1,i,1,1)
 120       athet(2,i,-1,1)=-athet(2,i,1,1)
 121       bthet(1,i,-1,1)=bthet(1,i,1,1)
 122       bthet(2,i,-1,1)=bthet(2,i,1,1)
 123       enddo
 124       do i=-ntyp,-1
 125       a0thet(i)=a0thet(-i)
 126       athet(1,i,-1,-1)=athet(1,-i,1,1)
 127       athet(2,i,-1,-1)=-athet(2,-i,1,1)
 128       bthet(1,i,-1,-1)=bthet(1,-i,1,1)
 129       bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
 130       athet(1,i,-1,1)=athet(1,-i,1,1)
 131       athet(2,i,-1,1)=-athet(2,-i,1,1)
 132       bthet(1,i,-1,1)=-bthet(1,-i,1,1)
 133       bthet(2,i,-1,1)=bthet(2,-i,1,1)
 134       athet(1,i,1,-1)=-athet(1,-i,1,1)
 135       athet(2,i,1,-1)=athet(2,-i,1,1)
 136       bthet(1,i,1,-1)=bthet(1,-i,1,1)
 137       bthet(2,i,1,-1)=-bthet(2,-i,1,1)
 138       theta0(i)=theta0(-i)
 139       sig0(i)=sig0(-i)
 140       sigc0(i)=sigc0(-i)
 141        do j=0,3
 142         polthet(j,i)=polthet(j,-i)
 143        enddo
 144        do j=1,3
 145          gthet(j,i)=gthet(j,-i)
 146        enddo
 147       enddo
 148
 149       close (ithep)
 150       if (lprint) then
 151       if (.not.LaTeX) then
 152         write (iout,'(a)')
 153      &    'Parameters of the virtual-bond valence angles:'
 154         write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
 155      & '    ATHETA0   ','         A1   ','        A2    ',
 156      & '        B1    ','         B2   '
 157         do i=1,ntyp
 158           write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
 159      &        a0thet(i),(athet(j,i,1,1),j=1,2),(bthet(j,i,1,1),j=1,2)
 160         enddo
 161         write (iout,'(/a/9x,5a/79(1h-))')
 162      & 'Parameters of the expression for sigma(theta_c):',
 163      & '     ALPH0    ','      ALPH1   ','     ALPH2    ',
 164      & '     ALPH3    ','    SIGMA0C   '
 165         do i=1,ntyp
 166           write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
 167      &      (polthet(j,i),j=0,3),sigc0(i)
 168         enddo
 169         write (iout,'(/a/9x,5a/79(1h-))')
 170      & 'Parameters of the second gaussian:',
 171      & '    THETA0    ','     SIGMA0   ','        G1    ',
 172      & '        G2    ','         G3   '
 173         do i=1,ntyp
 174           write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
 175      &       sig0(i),(gthet(j,i),j=1,3)
 176         enddo
 177        else
 178         write (iout,'(a)')
 179      &    'Parameters of the virtual-bond valence angles:'
 180         write (iout,'(/a/9x,5a/79(1h-))')
 181      & 'Coefficients of expansion',
 182      & '     theta0   ','    a1*10^2   ','   a2*10^2    ',
 183      & '   b1*10^1    ','    b2*10^1   '
 184         do i=1,ntyp
 185           write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),
 186      &        a0thet(i),(100*athet(j,i,1,1),j=1,2),
 187      &        (10*bthet(j,i,1,1),j=1,2)
 188         enddo
 189         write (iout,'(/a/9x,5a/79(1h-))')
 190      & 'Parameters of the expression for sigma(theta_c):',
 191      & ' alpha0       ','  alph1       ',' alph2        ',
 192      & ' alhp3        ','   sigma0c    '
 193         do i=1,ntyp
 194           write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),
 195      &      (polthet(j,i),j=0,3),sigc0(i)
 196         enddo
 197         write (iout,'(/a/9x,5a/79(1h-))')
 198      & 'Parameters of the second gaussian:',
 199      & '    theta0    ','  sigma0*10^2 ','      G1*10^-1',
 200      & '        G2    ','   G3*10^1    '
 201         do i=1,ntyp
 202           write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),
 203      &       100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
 204         enddo
 205       endif
 206       endif
 207 #else
 208 C
 209 C Read the parameters of Utheta determined from ab initio surfaces
 210 C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
 211 C
 212       read (ithep,*,err=111,end=111) nthetyp,ntheterm,ntheterm2,
 213      &  ntheterm3,nsingle,ndouble
 214       nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
 215       read (ithep,*,err=111,end=111) (ithetyp(i),i=1,ntyp1)
 216       do i=-ntyp1,-1
 217         ithetyp(i)=-ithetyp(-i)
 218       enddo
 219       do iblock=1,2
 220       do i=-maxthetyp,maxthetyp
 221         do j=-maxthetyp,maxthetyp
 222           do k=-maxthetyp,maxthetyp
 223             aa0thet(i,j,k,iblock)=0.0d0
 224             do l=1,ntheterm
 225               aathet(l,i,j,k,iblock)=0.0d0
 226             enddo
 227             do l=1,ntheterm2
 228               do m=1,nsingle
 229                 bbthet(m,l,i,j,k,iblock)=0.0d0
 230                 ccthet(m,l,i,j,k,iblock)=0.0d0
 231                 ddthet(m,l,i,j,k,iblock)=0.0d0
 232                 eethet(m,l,i,j,k,iblock)=0.0d0
 233               enddo
 234             enddo
 235             do l=1,ntheterm3
 236               do m=1,ndouble
 237                 do mm=1,ndouble
 238                  ffthet(mm,m,l,i,j,k,iblock)=0.0d0
 239                  ggthet(mm,m,l,i,j,k,iblock)=0.0d0
 240                 enddo
 241               enddo
 242             enddo
 243           enddo
 244         enddo
 245       enddo
 246       enddo
 247 c VAR:iblock means terminally blocking group 1=non-proline 2=proline
 248       do iblock=1,2
 249 c VAR:ntethtyp is type of theta potentials type currently 0=glycine
 250 c VAR:1=non-glicyne non-proline 2=proline
 251 c VAR:negative values for D-aminoacid
 252       do i=0,nthetyp
 253         do j=-nthetyp,nthetyp
 254           do k=-nthetyp,nthetyp
 255             read (ithep,'(6a)',end=111,err=111) res1
 256             read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock)
 257 c VAR: aa0thet is variable describing the average value of Foureir
 258 c VAR: expansion series
 259 c VAR: aathet is foureir expansion in theta/2 angle for full formula
 260 c VAR: look at the fitting equation in Kozlowska et al., J. Phys.:
 261 Condens. Matter 19 (2007) 285203 and Sieradzan et al., unpublished
 262             read (ithep,*,end=111,err=111)
 263      &(aathet(l,i,j,k,iblock),l=1,ntheterm)
 264             read (ithep,*,end=111,err=111)
 265      &       ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
 266      &        (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
 267      &        (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
 268      &        (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle),
 269      &        ll=1,ntheterm2)
 270             read (ithep,*,end=111,err=111)
 271      &      (((ffthet(llll,lll,ll,i,j,k,iblock),
 272      &         ffthet(lll,llll,ll,i,j,k,iblock),
 273      &         ggthet(llll,lll,ll,i,j,k,iblock),
 274      &         ggthet(lll,llll,ll,i,j,k,iblock),
 275      &         llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
 276           enddo
 277         enddo
 278       enddo
 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",ithep_pdb
 393       do i=1,ntyp
 394 c      write (2,*) 'i=',i
 395         read (ithep_pdb,*,err=111,end=111)
 396      &     a0thet(i),(athet(j,i,1,1),j=1,2),
 397      &    (bthet(j,i,1,1),j=1,2)
 398         read (ithep_pdb,*,err=111,end=111) (polthet(j,i),j=0,3)
 399         read (ithep_pdb,*,err=111,end=111) (gthet(j,i),j=1,3)
 400         read (ithep_pdb,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
 401         sigc0(i)=sigc0(i)**2
 402       enddo
 403       do i=1,ntyp
 404       athet(1,i,1,-1)=athet(1,i,1,1)
 405       athet(2,i,1,-1)=athet(2,i,1,1)
 406       bthet(1,i,1,-1)=-bthet(1,i,1,1)
 407       bthet(2,i,1,-1)=-bthet(2,i,1,1)
 408       athet(1,i,-1,1)=-athet(1,i,1,1)
 409       athet(2,i,-1,1)=-athet(2,i,1,1)
 410       bthet(1,i,-1,1)=bthet(1,i,1,1)
 411       bthet(2,i,-1,1)=bthet(2,i,1,1)
 412       enddo
 413       do i=-ntyp,-1
 414       a0thet(i)=a0thet(-i)
 415       athet(1,i,-1,-1)=athet(1,-i,1,1)
 416       athet(2,i,-1,-1)=-athet(2,-i,1,1)
 417       bthet(1,i,-1,-1)=bthet(1,-i,1,1)
 418       bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
 419       athet(1,i,-1,1)=athet(1,-i,1,1)
 420       athet(2,i,-1,1)=-athet(2,-i,1,1)
 421       bthet(1,i,-1,1)=-bthet(1,-i,1,1)
 422       bthet(2,i,-1,1)=bthet(2,-i,1,1)
 423       athet(1,i,1,-1)=-athet(1,-i,1,1)
 424       athet(2,i,1,-1)=athet(2,-i,1,1)
 425       bthet(1,i,1,-1)=bthet(1,-i,1,1)
 426       bthet(2,i,1,-1)=-bthet(2,-i,1,1)
 427       theta0(i)=theta0(-i)
 428       sig0(i)=sig0(-i)
 429       sigc0(i)=sigc0(-i)
 430        do j=0,3
 431         polthet(j,i)=polthet(j,-i)
 432        enddo
 433        do j=1,3
 434          gthet(j,i)=gthet(j,-i)
 435        enddo
 436       enddo
 437       write (2,*) "End reading THETA_PDB"
 438       close (ithep_pdb)
 439 #endif
 440       close(ithep)
 441 #ifdef CRYST_SC
 442 C
 443 C Read the parameters of the probability distribution/energy expression
 444 C of the side chains.
 445 C
 446       do i=1,ntyp
 447         read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
 448         if (i.eq.10) then
 449           dsc_inv(i)=0.0D0
 450         else
 451           dsc_inv(i)=1.0D0/dsc(i)
 452         endif
 453         if (i.ne.10) then
 454         do j=1,nlob(i)
 455           do k=1,3
 456             do l=1,3
 457               blower(l,k,j)=0.0D0
 458             enddo
 459           enddo
 460         enddo
 461         bsc(1,i)=0.0D0
 462         read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),
 463      &    ((blower(k,l,1),l=1,k),k=1,3)
 464         censc(1,1,-i)=censc(1,1,i)
 465         censc(2,1,-i)=censc(2,1,i)
 466         censc(3,1,-i)=-censc(3,1,i)
 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       do i=1,ntyp
 547         read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
 548         if (i.eq.10) then
 549           dsc_inv(i)=0.0D0
 550         else
 551           dsc_inv(i)=1.0D0/dsc(i)
 552         endif
 553         if (i.ne.10) then
 554         do j=1,nlob(i)
 555           do k=1,3
 556             do l=1,3
 557               blower(l,k,j)=0.0D0
 558             enddo
 559           enddo
 560         enddo
 561         bsc(1,i)=0.0D0
 562         read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),
 563      &    ((blower(k,l,1),l=1,k),k=1,3)
 564         do j=2,nlob(i)
 565           read (irotam_pdb,*,end=112,err=112) bsc(j,i)
 566           read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),
 567      &                                 ((blower(k,l,j),l=1,k),k=1,3)
 568         enddo
 569         do j=1,nlob(i)
 570           do k=1,3
 571             do l=1,k
 572               akl=0.0D0
 573               do m=1,3
 574                 akl=akl+blower(k,m,j)*blower(l,m,j)
 575               enddo
 576               gaussc(k,l,j,i)=akl
 577               gaussc(l,k,j,i)=akl
 578             enddo
 579           enddo
 580         enddo
 581         endif
 582       enddo
 583       close (irotam_pdb)
 584       write (2,*) "End reading ROTAM_PDB"
 585 #endif
 586       close(irotam)
 587
 588 #ifdef CRYST_TOR
 589 C
 590 C Read torsional parameters in old format
 591 C
 592       read (itorp,*,end=113,err=113) ntortyp,nterm_old
 593       if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
 594       read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
 595       do i=1,ntortyp
 596         do j=1,ntortyp
 597           read (itorp,'(a)')
 598           do k=1,nterm_old
 599             read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
 600           enddo
 601         enddo
 602       enddo
 603       close (itorp)
 604       if (lprint) then
 605         write (iout,'(/a/)') 'Torsional constants:'
 606         do i=1,ntortyp
 607           do j=1,ntortyp
 608             write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
 609             write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
 610           enddo
 611         enddo
 612       endif
 613 #else
 614 C
 615 C Read torsional parameters
 616 C
 617       read (itorp,*,end=113,err=113) ntortyp
 618       read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
 619       do iblock=1,2
 620       do i=-ntyp,-1
 621        itortyp(i)=-itortyp(-i)
 622       enddo
 623       write (iout,*) 'ntortyp',ntortyp
 624       do i=0,ntortyp-1
 625         do j=-ntortyp+1,ntortyp-1
 626           read (itorp,*,end=113,err=113) nterm(i,j,iblock),
 627      &          nlor(i,j,iblock)
 628           nterm(-i,-j,iblock)=nterm(i,j,iblock)
 629           nlor(-i,-j,iblock)=nlor(i,j,iblock)
 630           v0ij=0.0d0
 631           si=-1.0d0
 632           do k=1,nterm(i,j,iblock)
 633             read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),
 634      &      v2(k,i,j,iblock)
 635             v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
 636             v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
 637             v0ij=v0ij+si*v1(k,i,j,iblock)
 638             si=-si
 639 c         write(iout,*) i,j,k,iblock,nterm(i,j,iblock)
 640 c         write(iout,*) v1(k,-i,-j,iblock),v1(k,i,j,iblock),
 641 c      &v2(k,-i,-j,iblock),v2(k,i,j,iblock)
 642           enddo
 643           do k=1,nlor(i,j,iblock)
 644             read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),
 645      &        vlor2(k,i,j),vlor3(k,i,j)
 646             v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
 647           enddo
 648           v0(i,j,iblock)=v0ij
 649           v0(-i,-j,iblock)=v0ij
 650         enddo
 651       enddo
 652       enddo
 653       close (itorp)
 654       if (lprint) then
 655         write (iout,'(/a/)') 'Torsional constants:'
 656         do i=1,ntortyp
 657           do j=1,ntortyp
 658             write (iout,*) 'ityp',i,' jtyp',j
 659             write (iout,*) 'Fourier constants'
 660             do k=1,nterm(i,j,iblock)
 661               write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
 662      &        v2(k,i,j,iblock)
 663             enddo
 664             write (iout,*) 'Lorenz constants'
 665             do k=1,nlor(i,j,iblock)
 666               write (iout,'(3(1pe15.5))')
 667      &         vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
 668             enddo
 669           enddo
 670         enddo
 671       endif
 672
 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)
 788      &nterm_sccor(i,j),nlor_sccor(i,j)
 789           v0ijsccor=0.0d0
 790           v0ijsccor1=0.0d0
 791           v0ijsccor2=0.0d0
 792           v0ijsccor3=0.0d0
 793           si=-1.0d0
 794           nterm_sccor(-i,j)=nterm_sccor(i,j)
 795           nterm_sccor(-i,-j)=nterm_sccor(i,j)
 796           nterm_sccor(i,-j)=nterm_sccor(i,j)
 797           do k=1,nterm_sccor(i,j)
 798             read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
 799      &    ,v2sccor(k,l,i,j)
 800             if (j.eq.iscprol) then
 801              if (i.eq.isccortyp(10)) then
 802              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
 803              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
 804              else
 805              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
 806      &                        +v2sccor(k,l,i,j)*dsqrt(0.75d0)
 807              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
 808      &                        +v1sccor(k,l,i,j)*dsqrt(0.75d0)
 809              v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
 810              v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
 811              v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
 812              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
 813              endif
 814             else
 815              if (i.eq.isccortyp(10)) then
 816              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
 817              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
 818              else
 819                if (j.eq.isccortyp(10)) then
 820              v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
 821              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
 822                else
 823              v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
 824              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
 825              v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
 826              v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
 827              v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
 828              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
 829                 endif
 830                endif
 831             endif
 832             v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
 833             v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
 834             v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
 835             v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
 836             si=-si
 837           enddo
 838           do k=1,nlor_sccor(i,j)
 839             read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
 840      &        vlor2sccor(k,i,j),vlor3sccor(k,i,j)
 841             v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
 842      &(1+vlor3sccor(k,i,j)**2)
 843           enddo
 844           v0sccor(l,i,j)=v0ijsccor
 845           v0sccor(l,-i,j)=v0ijsccor1
 846           v0sccor(l,i,-j)=v0ijsccor2
 847           v0sccor(l,-i,-j)=v0ijsccor3
 848         enddo
 849       enddo
 850       enddo
 851       close (isccor)
 852 #else
 853       read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
 854 c      write (iout,*) 'ntortyp',ntortyp
 855       maxinter=3
 856 cc maxinter is maximum interaction sites
 857       do l=1,maxinter
 858       do i=1,nsccortyp
 859         do j=1,nsccortyp
 860           read (isccor,*,end=119,err=119)
 861      & nterm_sccor(i,j),nlor_sccor(i,j)
 862           v0ijsccor=0.0d0
 863           si=-1.0d0
 864
 865           do k=1,nterm_sccor(i,j)
 866             read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
 867      &    ,v2sccor(k,l,i,j)
 868             v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
 869             si=-si
 870           enddo
 871           do k=1,nlor_sccor(i,j)
 872             read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
 873      &        vlor2sccor(k,i,j),vlor3sccor(k,i,j)
 874             v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
 875      &(1+vlor3sccor(k,i,j)**2)
 876           enddo
 877           v0sccor(i,j,iblock)=v0ijsccor
 878         enddo
 879       enddo
 880       enddo
 881       close (isccor)
 882
 883 #endif
 884       if (lprint) then
 885         write (iout,'(/a/)') 'Torsional constants:'
 886         do i=1,nsccortyp
 887           do j=1,nsccortyp
 888             write (iout,*) 'ityp',i,' jtyp',j
 889             write (iout,*) 'Fourier constants'
 890             do k=1,nterm_sccor(i,j)
 891       write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
 892             enddo
 893             write (iout,*) 'Lorenz constants'
 894             do k=1,nlor_sccor(i,j)
 895               write (iout,'(3(1pe15.5))')
 896      &         vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
 897             enddo
 898           enddo
 899         enddo
 900       endif
 901
 902 C
 903 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
 904 C         interaction energy of the Gly, Ala, and Pro prototypes.
 905 C
 906       if (lprint) then
 907         write (iout,*)
 908         write (iout,*) "Coefficients of the cumulants"
 909       endif
 910       read (ifourier,*) nloctyp
 911       do i=0,nloctyp-1
 912         read (ifourier,*,end=115,err=115)
 913         read (ifourier,*,end=115,err=115) (b(ii),ii=1,13)
 914         if (lprint) then
 915         write (iout,*) 'Type',i
 916         write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
 917         endif
 918         B1(1,i)  = b(3)
 919         B1(2,i)  = b(5)
 920         B1(1,-i) = b(3)
 921         B1(2,-i) = -b(5)
 922 c        b1(1,i)=0.0d0
 923 c        b1(2,i)=0.0d0
 924         B1tilde(1,i) = b(3)
 925         B1tilde(2,i) =-b(5)
 926         B1tilde(1,-i) =-b(3)
 927         B1tilde(2,-i) =b(5)
 928         b1tilde(1,i)=0.0d0
 929         b1tilde(2,i)=0.0d0
 930         B2(1,i)  = b(2)
 931         B2(2,i)  = b(4)
 932         B2(1,-i)  =b(2)
 933         B2(2,-i)  =-b(4)
 934
 935 c        b2(1,i)=0.0d0
 936 c        b2(2,i)=0.0d0
 937         CC(1,1,i)= b(7)
 938         CC(2,2,i)=-b(7)
 939         CC(2,1,i)= b(9)
 940         CC(1,2,i)= b(9)
 941         CC(1,1,-i)= b(7)
 942         CC(2,2,-i)=-b(7)
 943         CC(2,1,-i)=-b(9)
 944         CC(1,2,-i)=-b(9)
 945 c        CC(1,1,i)=0.0d0
 946 c        CC(2,2,i)=0.0d0
 947 c        CC(2,1,i)=0.0d0
 948 c        CC(1,2,i)=0.0d0
 949         Ctilde(1,1,i)=b(7)
 950         Ctilde(1,2,i)=b(9)
 951         Ctilde(2,1,i)=-b(9)
 952         Ctilde(2,2,i)=b(7)
 953         Ctilde(1,1,-i)=b(7)
 954         Ctilde(1,2,-i)=-b(9)
 955         Ctilde(2,1,-i)=b(9)
 956         Ctilde(2,2,-i)=b(7)
 957
 958 c        Ctilde(1,1,i)=0.0d0
 959 c        Ctilde(1,2,i)=0.0d0
 960 c        Ctilde(2,1,i)=0.0d0
 961 c        Ctilde(2,2,i)=0.0d0
 962         DD(1,1,i)= b(6)
 963         DD(2,2,i)=-b(6)
 964         DD(2,1,i)= b(8)
 965         DD(1,2,i)= b(8)
 966         DD(1,1,-i)= b(6)
 967         DD(2,2,-i)=-b(6)
 968         DD(2,1,-i)=-b(8)
 969         DD(1,2,-i)=-b(8)
 970 c        DD(1,1,i)=0.0d0
 971 c        DD(2,2,i)=0.0d0
 972 c        DD(2,1,i)=0.0d0
 973 c        DD(1,2,i)=0.0d0
 974         Dtilde(1,1,i)=b(6)
 975         Dtilde(1,2,i)=b(8)
 976         Dtilde(2,1,i)=-b(8)
 977         Dtilde(2,2,i)=b(6)
 978         Dtilde(1,1,-i)=b(6)
 979         Dtilde(1,2,-i)=-b(8)
 980         Dtilde(2,1,-i)=b(8)
 981         Dtilde(2,2,-i)=b(6)
 982
 983 c        Dtilde(1,1,i)=0.0d0
 984 c        Dtilde(1,2,i)=0.0d0
 985 c        Dtilde(2,1,i)=0.0d0
 986 c        Dtilde(2,2,i)=0.0d0
 987         EE(1,1,i)= b(10)+b(11)
 988         EE(2,2,i)=-b(10)+b(11)
 989         EE(2,1,i)= b(12)-b(13)
 990         EE(1,2,i)= b(12)+b(13)
 991         EE(1,1,-i)= b(10)+b(11)
 992         EE(2,2,-i)=-b(10)+b(11)
 993         EE(2,1,-i)=-b(12)+b(13)
 994         EE(1,2,-i)=-b(12)-b(13)
 995
 996 c        ee(1,1,i)=1.0d0
 997 c        ee(2,2,i)=1.0d0
 998 c        ee(2,1,i)=0.0d0
 999 c        ee(1,2,i)=0.0d0
1000 c        ee(2,1,i)=ee(1,2,i)
1001       enddo
1002       if (lprint) then
1003       do i=1,nloctyp
1004         write (iout,*) 'Type',i
1005         write (iout,*) 'B1'
1006         write(iout,*) B1(1,i),B1(2,i)
1007         write (iout,*) 'B2'
1008         write(iout,*) B2(1,i),B2(2,i)
1009         write (iout,*) 'CC'
1010         do j=1,2
1011           write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
1012         enddo
1013         write(iout,*) 'DD'
1014         do j=1,2
1015           write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
1016         enddo
1017         write(iout,*) 'EE'
1018         do j=1,2
1019           write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
1020         enddo
1021       enddo
1022       endif
1023
1024 C
1025 C Read electrostatic-interaction parameters
1026 C
1027       if (lprint) then
1028         write (iout,*)
1029         write (iout,'(/a)') 'Electrostatic interaction constants:'
1030         write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
1031      &            'IT','JT','APP','BPP','AEL6','AEL3'
1032       endif
1033       read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
1034       read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
1035       read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
1036       read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
1037       close (ielep)
1038       do i=1,2
1039         do j=1,2
1040         rri=rpp(i,j)**6
1041         app (i,j)=epp(i,j)*rri*rri
1042         bpp (i,j)=-2.0D0*epp(i,j)*rri
1043         ael6(i,j)=elpp6(i,j)*4.2D0**6
1044         ael3(i,j)=elpp3(i,j)*4.2D0**3
1045 c        lprint=.true.
1046         if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
1047      &                    ael6(i,j),ael3(i,j)
1048 c        lprint=.false.
1049         enddo
1050       enddo
1051 C
1052 C Read side-chain interaction parameters.
1053 C
1054       read (isidep,*,end=117,err=117) ipot,expon
1055       if (ipot.lt.1 .or. ipot.gt.5) then
1056         write (iout,'(2a)') 'Error while reading SC interaction',
1057      &               'potential file - unknown potential type.'
1058 #ifdef MPI
1059         call MPI_Finalize(Ierror)
1060 #endif
1061         stop
1062       endif
1063       expon2=expon/2
1064       if(me.eq.king)
1065      & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
1066      & ', exponents are ',expon,2*expon
1067       goto (10,20,30,30,40) ipot
1068 C----------------------- LJ potential ---------------------------------
1069    10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1070      &   (sigma0(i),i=1,ntyp)
1071       if (lprint) then
1072         write (iout,'(/a/)') 'Parameters of the LJ potential:'
1073         write (iout,'(a/)') 'The epsilon array:'
1074         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1075         write (iout,'(/a)') 'One-body parameters:'
1076         write (iout,'(a,4x,a)') 'residue','sigma'
1077         write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1078       endif
1079       goto 50
1080 C----------------------- LJK potential --------------------------------
1081    20 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1082      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1083       if (lprint) then
1084         write (iout,'(/a/)') 'Parameters of the LJK potential:'
1085         write (iout,'(a/)') 'The epsilon array:'
1086         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1087         write (iout,'(/a)') 'One-body parameters:'
1088         write (iout,'(a,4x,2a)') 'residue','   sigma  ','    r0    '
1089         write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
1090      &        i=1,ntyp)
1091       endif
1092       goto 50
1093 C---------------------- GB or BP potential -----------------------------
1094    30 do i=1,ntyp
1095        read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
1096       enddo
1097       read (isidep,*,end=117,err=117)(sigma0(i),i=1,ntyp)
1098       read (isidep,*,end=117,err=117)(sigii(i),i=1,ntyp)
1099       read (isidep,*,end=117,err=117)(chip(i),i=1,ntyp)
1100       read (isidep,*,end=117,err=117)(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=117,err=117)((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 c      lprint=.true.
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 c      lprint=.false.
1256 #endif
1257 C
1258 C Define the constants of the disulfide bridge
1259 C
1260 C      ebr=-12.00D0
1261 c
1262 c Old arbitrary potential - commented out.
1263 c
1264 c      dbr= 4.20D0
1265 c      fbr= 3.30D0
1266 c
1267 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1268 c energy surface of diethyl disulfide.
1269 c A. Liwo and U. Kozlowska, 11/24/03
1270 c
1271 C      D0CM = 3.78d0
1272 C      AKCM = 15.1d0
1273 C      AKTH = 11.0d0
1274 C      AKCT = 12.0d0
1275 C      V1SS =-1.08d0
1276 C      V2SS = 7.61d0
1277 C      V3SS = 13.7d0
1278 c      akcm=0.0d0
1279 c      akth=0.0d0
1280 c      akct=0.0d0
1281 c      v1ss=0.0d0
1282 c      v2ss=0.0d0
1283 c      v3ss=0.0d0
1284
1285 C      if(me.eq.king) then
1286 C      write (iout,'(/a)') "Disulfide bridge parameters:"
1287 C      write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1288 C      write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1289 C      write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1290 C      write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1291 C     &  ' v3ss:',v3ss
1292 C      endif
1293       return
1294   111 write (iout,*) "Error reading bending energy parameters."
1295       goto 999
1296   112 write (iout,*) "Error reading rotamer energy parameters."
1297       goto 999
1298   113 write (iout,*) "Error reading torsional energy parameters."
1299       goto 999
1300   114 write (iout,*) "Error reading double torsional energy parameters."
1301       goto 999
1302   115 write (iout,*)
1303      &  "Error reading cumulant (multibody energy) parameters."
1304       goto 999
1305   116 write (iout,*) "Error reading electrostatic energy parameters."
1306       goto 999
1307   117 write (iout,*) "Error reading side chain interaction parameters."
1308       goto 999
1309   118 write (iout,*) "Error reading SCp interaction parameters."
1310       goto 999
1311   119 write (iout,*) "Error reading SCCOR parameters"
1312   999 continue
1313 #ifdef MPI
1314       call MPI_Finalize(Ierror)
1315 #endif
1316       stop
1317       return
1318       end
1319
1320
1321       subroutine getenv_loc(var, val)
1322       character(*) var, val
1323
1324 #ifdef WINIFL
1325       character(2000) line
1326       external ilen
1327
1328       open (196,file='env',status='old',readonly,shared)
1329       iread=0
1330 c      write(*,*)'looking for ',var
1331 10    read(196,*,err=11,end=11)line
1332       iread=index(line,var)
1333 c      write(*,*)iread,' ',var,' ',line
1334       if (iread.eq.0) go to 10
1335 c      write(*,*)'---> ',line
1336 11    continue
1337       if(iread.eq.0) then
1338 c       write(*,*)'CHUJ'
1339        val=''
1340       else
1341        iread=iread+ilen(var)+1
1342        read (line(iread:),*,err=12,end=12) val
1343 c       write(*,*)'OK: ',var,' = ',val
1344       endif
1345       close(196)
1346       return
1347 12    val=''
1348       close(196)
1349 #elif (defined CRAY)
1350       integer lennam,lenval,ierror
1351 c
1352 c        getenv using a POSIX call, useful on the T3D
1353 c        Sept 1996, comment out error check on advice of H. Pritchard
1354 c
1355       lennam = len(var)
1356       if(lennam.le.0) stop '--error calling getenv--'
1357       call pxfgetenv(var,lennam,val,lenval,ierror)
1358 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1359 #else
1360       call getenv(var,val)
1361 #endif
1362
1363       return
1364       end