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,*,err=111,end=111) a0thet(i),(athet(j,i,1,1),j=1,2),
 395      &    (bthet(j,i,1,1),j=1,2)
 396         read (ithep,*,err=111,end=111) (polthet(j,i),j=0,3)
 397         read (ithep,*,err=111,end=111) (gthet(j,i),j=1,3)
 398         read (ithep,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
 399         sigc0(i)=sigc0(i)**2
 400       enddo
 401       do i=1,ntyp
 402       athet(1,i,1,-1)=athet(1,i,1,1)
 403       athet(2,i,1,-1)=athet(2,i,1,1)
 404       bthet(1,i,1,-1)=-bthet(1,i,1,1)
 405       bthet(2,i,1,-1)=-bthet(2,i,1,1)
 406       athet(1,i,-1,1)=-athet(1,i,1,1)
 407       athet(2,i,-1,1)=-athet(2,i,1,1)
 408       bthet(1,i,-1,1)=bthet(1,i,1,1)
 409       bthet(2,i,-1,1)=bthet(2,i,1,1)
 410       enddo
 411       do i=-ntyp,-1
 412       a0thet(i)=a0thet(-i)
 413       athet(1,i,-1,-1)=athet(1,-i,1,1)
 414       athet(2,i,-1,-1)=-athet(2,-i,1,1)
 415       bthet(1,i,-1,-1)=bthet(1,-i,1,1)
 416       bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
 417       athet(1,i,-1,1)=athet(1,-i,1,1)
 418       athet(2,i,-1,1)=-athet(2,-i,1,1)
 419       bthet(1,i,-1,1)=-bthet(1,-i,1,1)
 420       bthet(2,i,-1,1)=bthet(2,-i,1,1)
 421       athet(1,i,1,-1)=-athet(1,-i,1,1)
 422       athet(2,i,1,-1)=athet(2,-i,1,1)
 423       bthet(1,i,1,-1)=bthet(1,-i,1,1)
 424       bthet(2,i,1,-1)=-bthet(2,-i,1,1)
 425       theta0(i)=theta0(-i)
 426       sig0(i)=sig0(-i)
 427       sigc0(i)=sigc0(-i)
 428        do j=0,3
 429         polthet(j,i)=polthet(j,-i)
 430        enddo
 431        do j=1,3
 432          gthet(j,i)=gthet(j,-i)
 433        enddo
 434       enddo
 435       write (2,*) "End reading THETA_PDB"
 436       close (ithep_pdb)
 437 #endif
 438       close(ithep)
 439 #ifdef CRYST_SC
 440 C
 441 C Read the parameters of the probability distribution/energy expression
 442 C of the side chains.
 443 C
 444       do i=1,ntyp
 445         read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
 446         if (i.eq.10) then
 447           dsc_inv(i)=0.0D0
 448         else
 449           dsc_inv(i)=1.0D0/dsc(i)
 450         endif
 451         if (i.ne.10) then
 452         do j=1,nlob(i)
 453           do k=1,3
 454             do l=1,3
 455               blower(l,k,j)=0.0D0
 456             enddo
 457           enddo
 458         enddo
 459         bsc(1,i)=0.0D0
 460         read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),
 461      &    ((blower(k,l,1),l=1,k),k=1,3)
 462         censc(1,1,-i)=censc(1,1,i)
 463         censc(2,1,-i)=censc(2,1,i)
 464         censc(3,1,-i)=-censc(3,1,i)
 465
 466         do j=2,nlob(i)
 467           read (irotam,*,end=112,err=112) bsc(j,i)
 468           read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),
 469      &                                 ((blower(k,l,j),l=1,k),k=1,3)
 470         censc(1,j,-i)=censc(1,j,i)
 471         censc(2,j,-i)=censc(2,j,i)
 472         censc(3,j,-i)=-censc(3,j,i)
 473 C BSC is amplitude of Gaussian
 474         enddo
 475         do j=1,nlob(i)
 476           do k=1,3
 477             do l=1,k
 478               akl=0.0D0
 479               do m=1,3
 480                 akl=akl+blower(k,m,j)*blower(l,m,j)
 481               enddo
 482               gaussc(k,l,j,i)=akl
 483               gaussc(l,k,j,i)=akl
 484               if (((k.eq.3).and.(l.ne.3))
 485      &        .or.((l.eq.3).and.(k.ne.3))) then
 486                 gaussc(k,l,j,-i)=-akl
 487                 gaussc(l,k,j,-i)=-akl
 488               else
 489                 gaussc(k,l,j,-i)=akl
 490                 gaussc(l,k,j,-i)=akl
 491               endif
 492             enddo
 493           enddo
 494         enddo
 495         endif
 496       enddo
 497       close (irotam)
 498       if (lprint) then
 499         write (iout,'(/a)') 'Parameters of side-chain local geometry'
 500         do i=1,ntyp
 501           nlobi=nlob(i)
 502           if (nlobi.gt.0) then
 503             if (LaTeX) then
 504               write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
 505      &         ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
 506                write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
 507      &                             'log h',(bsc(j,i),j=1,nlobi)
 508                write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
 509      &        'x',((censc(k,j,i),k=1,3),j=1,nlobi)
 510               do k=1,3
 511                 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
 512      &                 ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
 513               enddo
 514             else
 515               write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
 516               write (iout,'(a,f10.4,4(16x,f10.4))')
 517      &                             'Center  ',(bsc(j,i),j=1,nlobi)
 518               write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),
 519      &           j=1,nlobi)
 520               write (iout,'(a)')
 521             endif
 522           endif
 523         enddo
 524       endif
 525 #else
 526 C
 527 C Read scrot parameters for potentials determined from all-atom AM1 calculations
 528 C added by Urszula Kozlowska 07/11/2007
 529 C
 530       do i=1,ntyp
 531         read (irotam,*,end=112,err=112)
 532        if (i.eq.10) then
 533          read (irotam,*,end=112,err=112)
 534        else
 535          do j=1,65
 536            read(irotam,*,end=112,err=112) sc_parmin(j,i)
 537          enddo
 538        endif
 539       enddo
 540 C
 541 C Read the parameters of the probability distribution/energy expression
 542 C of the side chains.
 543 C
 544       do i=1,ntyp
 545         read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
 546         if (i.eq.10) then
 547           dsc_inv(i)=0.0D0
 548         else
 549           dsc_inv(i)=1.0D0/dsc(i)
 550         endif
 551         if (i.ne.10) then
 552         do j=1,nlob(i)
 553           do k=1,3
 554             do l=1,3
 555               blower(l,k,j)=0.0D0
 556             enddo
 557           enddo
 558         enddo
 559         bsc(1,i)=0.0D0
 560         read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),
 561      &    ((blower(k,l,1),l=1,k),k=1,3)
 562         do j=2,nlob(i)
 563           read (irotam_pdb,*,end=112,err=112) bsc(j,i)
 564           read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),
 565      &                                 ((blower(k,l,j),l=1,k),k=1,3)
 566         enddo
 567         do j=1,nlob(i)
 568           do k=1,3
 569             do l=1,k
 570               akl=0.0D0
 571               do m=1,3
 572                 akl=akl+blower(k,m,j)*blower(l,m,j)
 573               enddo
 574               gaussc(k,l,j,i)=akl
 575               gaussc(l,k,j,i)=akl
 576             enddo
 577           enddo
 578         enddo
 579         endif
 580       enddo
 581       close (irotam_pdb)
 582 #endif
 583       close(irotam)
 584
 585 #ifdef CRYST_TOR
 586 C
 587 C Read torsional parameters in old format
 588 C
 589       read (itorp,*,end=113,err=113) ntortyp,nterm_old
 590       if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
 591       read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
 592       do i=1,ntortyp
 593         do j=1,ntortyp
 594           read (itorp,'(a)')
 595           do k=1,nterm_old
 596             read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
 597           enddo
 598         enddo
 599       enddo
 600       close (itorp)
 601       if (lprint) then
 602         write (iout,'(/a/)') 'Torsional constants:'
 603         do i=1,ntortyp
 604           do j=1,ntortyp
 605             write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
 606             write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
 607           enddo
 608         enddo
 609       endif
 610 #else
 611 C
 612 C Read torsional parameters
 613 C
 614       read (itorp,*,end=113,err=113) ntortyp
 615       read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
 616       do iblock=1,2
 617       do i=-ntyp,-1
 618        itortyp(i)=-itortyp(-i)
 619       enddo
 620 c      write (iout,*) 'ntortyp',ntortyp
 621       do i=0,ntortyp-1
 622         do j=-ntortyp+1,ntortyp-1
 623           read (itorp,*,end=113,err=113) nterm(i,j,iblock),
 624      &          nlor(i,j,iblock)
 625           nterm(-i,-j,iblock)=nterm(i,j,iblock)
 626           nlor(-i,-j,iblock)=nlor(i,j,iblock)
 627           v0ij=0.0d0
 628           si=-1.0d0
 629           do k=1,nterm(i,j,iblock)
 630             read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),
 631      &      v2(k,i,j,iblock)
 632             v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
 633             v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
 634             v0ij=v0ij+si*v1(k,i,j,iblock)
 635             si=-si
 636 c         write(iout,*) i,j,k,iblock,nterm(i,j,iblock)
 637 c         write(iout,*) v1(k,-i,-j,iblock),v1(k,i,j,iblock),
 638 c      &v2(k,-i,-j,iblock),v2(k,i,j,iblock)
 639           enddo
 640           do k=1,nlor(i,j,iblock)
 641             read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),
 642      &        vlor2(k,i,j),vlor3(k,i,j)
 643             v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
 644           enddo
 645           v0(i,j,iblock)=v0ij
 646           v0(-i,-j,iblock)=v0ij
 647         enddo
 648       enddo
 649       enddo
 650       close (itorp)
 651       if (lprint) then
 652         write (iout,'(/a/)') 'Torsional constants:'
 653         do i=1,ntortyp
 654           do j=1,ntortyp
 655             write (iout,*) 'ityp',i,' jtyp',j
 656             write (iout,*) 'Fourier constants'
 657             do k=1,nterm(i,j,iblock)
 658               write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
 659      &        v2(k,i,j,iblock)
 660             enddo
 661             write (iout,*) 'Lorenz constants'
 662             do k=1,nlor(i,j,iblock)
 663               write (iout,'(3(1pe15.5))')
 664      &         vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
 665             enddo
 666           enddo
 667         enddo
 668       endif
 669 C
 670 C 6/23/01 Read parameters for double torsionals
 671 C
 672       do iblock=1,2
 673       do i=0,ntortyp-1
 674         do j=-ntortyp+1,ntortyp-1
 675           do k=-ntortyp+1,ntortyp-1
 676             read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
 677 c              write (iout,*) "OK onelett",
 678 c     &         i,j,k,t1,t2,t3
 679
 680             if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
 681      &        .or. t3.ne.toronelet(k)) then
 682               write (iout,*) "Error in double torsional parameter file",
 683      &         i,j,k,t1,t2,t3
 684 #ifdef MPI
 685               call MPI_Finalize(Ierror)
 686 #endif
 687                stop "Error in double torsional parameter file"
 688             endif
 689             read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),
 690      &         ntermd_2(i,j,k,iblock)
 691             ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
 692             ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
 693             read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,
 694      &         ntermd_1(i,j,k,iblock))
 695             read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,
 696      &         ntermd_1(i,j,k,iblock))
 697             read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,
 698      &         ntermd_1(i,j,k,iblock))
 699             read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,
 700      &         ntermd_1(i,j,k,iblock))
 701 C Martix of D parameters for one dimesional foureir series
 702             do l=1,ntermd_1(i,j,k,iblock)
 703              v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
 704              v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
 705              v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
 706              v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
 707 c            write(iout,*) "whcodze" ,
 708 c     & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
 709             enddo
 710             read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),
 711      &         v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
 712      &         v2s(m,l,i,j,k,iblock),
 713      &         m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
 714 C Martix of D parameters for two dimesional fourier series
 715             do l=1,ntermd_2(i,j,k,iblock)
 716              do m=1,l-1
 717              v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
 718              v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
 719              v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
 720              v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
 721              enddo!m
 722             enddo!l
 723           enddo!k
 724         enddo!j
 725       enddo!i
 726       enddo!iblock
 727       if (lprint) then
 728       write (iout,*)
 729       write (iout,*) 'Constants for double torsionals'
 730       do iblock=1,2
 731       do i=0,ntortyp-1
 732         do j=-ntortyp+1,ntortyp-1
 733           do k=-ntortyp+1,ntortyp-1
 734             write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
 735      &        ' nsingle',ntermd_1(i,j,k,iblock),
 736      &        ' ndouble',ntermd_2(i,j,k,iblock)
 737             write (iout,*)
 738             write (iout,*) 'Single angles:'
 739             do l=1,ntermd_1(i,j,k,iblock)
 740               write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
 741      &           v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
 742      &           v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),
 743      &           v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
 744             enddo
 745             write (iout,*)
 746             write (iout,*) 'Pairs of angles:'
 747             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
 748             do l=1,ntermd_2(i,j,k,iblock)
 749               write (iout,'(i5,20f10.5)')
 750      &         l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
 751             enddo
 752             write (iout,*)
 753             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
 754             do l=1,ntermd_2(i,j,k,iblock)
 755               write (iout,'(i5,20f10.5)')
 756      &         l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),
 757      &         (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
 758             enddo
 759             write (iout,*)
 760           enddo
 761         enddo
 762       enddo
 763       enddo
 764       endif
 765 #endif
 766 C Read of Side-chain backbone correlation parameters
 767 C Modified 11 May 2012 by Adasko
 768 CCC
 769 C
 770       read (isccor,*,end=119,err=119) nsccortyp
 771 #ifdef SCCORPDB
 772       read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
 773       do i=-ntyp,-1
 774         isccortyp(i)=-isccortyp(-i)
 775       enddo
 776       iscprol=isccortyp(20)
 777 c      write (iout,*) 'ntortyp',ntortyp
 778       maxinter=3
 779 cc maxinter is maximum interaction sites
 780       do l=1,maxinter
 781       do i=1,nsccortyp
 782         do j=1,nsccortyp
 783           read (isccor,*,end=119,err=119) nterm_sccor(i,j),nlor_sccor(i,j)
 784           v0ijsccor=0.0d0
 785           v0ijsccor1=0.0d0
 786           v0ijsccor2=0.0d0
 787           v0ijsccor3=0.0d0
 788           si=-1.0d0
 789           nterm_sccor(-i,j)=nterm_sccor(i,j)
 790           nterm_sccor(-i,-j)=nterm_sccor(i,j)
 791           nterm_sccor(i,-j)=nterm_sccor(i,j)
 792           do k=1,nterm_sccor(i,j)
 793             read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
 794      &    ,v2sccor(k,l,i,j)
 795             if (j.eq.iscprol) then
 796               if (i.eq.isccortyp(10)) then
 797               v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
 798               v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
 799               else
 800              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
 801      &                        +v2sccor(k,l,i,j)*dsqrt(0.75d0)
 802              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
 803      &                        +v1sccor(k,l,i,j)*dsqrt(0.75d0)
 804              v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
 805              v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
 806              v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
 807              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
 808              endif
 809             else
 810               if (i.eq.isccortyp(10)) then
 811               v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
 812               v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
 813               else
 814                 if (j.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              v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
 819              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
 820              v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
 821              v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
 822              v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
 823              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
 824             endif
 825              endif
 826              endif
 827             v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
 828             v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
 829             v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
 830             v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
 831             si=-si
 832           enddo
 833           do k=1,nlor_sccor(i,j)
 834             read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
 835      &        vlor2sccor(k,i,j),vlor3sccor(k,i,j)
 836             v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
 837      &(1+vlor3sccor(k,i,j)**2)
 838           enddo
 839           v0sccor(l,i,j)=v0ijsccor
 840           v0sccor(l,-i,j)=v0ijsccor1
 841           v0sccor(l,i,-j)=v0ijsccor2
 842           v0sccor(l,-i,-j)=v0ijsccor3
 843         enddo
 844       enddo
 845       enddo
 846       close (isccor)
 847 #else
 848       read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
 849 c      write (iout,*) 'ntortyp',ntortyp
 850       maxinter=3
 851 cc maxinter is maximum interaction sites
 852       do l=1,maxinter
 853       do i=1,nsccortyp
 854         do j=1,nsccortyp
 855           read (isccor,*,end=119,err=119)
 856      & nterm_sccor(i,j),nlor_sccor(i,j)
 857           v0ijsccor=0.0d0
 858           si=-1.0d0
 859
 860           do k=1,nterm_sccor(i,j)
 861             read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
 862      &    ,v2sccor(k,l,i,j)
 863             v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
 864             si=-si
 865           enddo
 866           do k=1,nlor_sccor(i,j)
 867             read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
 868      &        vlor2sccor(k,i,j),vlor3sccor(k,i,j)
 869             v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
 870      &(1+vlor3sccor(k,i,j)**2)
 871           enddo
 872           v0sccor(i,j,iblock)=v0ijsccor
 873         enddo
 874       enddo
 875       enddo
 876       close (isccor)
 877
 878 #endif
 879       if (lprint) then
 880         write (iout,'(/a/)') 'Torsional constants:'
 881         do i=1,nsccortyp
 882           do j=1,nsccortyp
 883             write (iout,*) 'ityp',i,' jtyp',j
 884             write (iout,*) 'Fourier constants'
 885             do k=1,nterm_sccor(i,j)
 886       write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
 887             enddo
 888             write (iout,*) 'Lorenz constants'
 889             do k=1,nlor_sccor(i,j)
 890               write (iout,'(3(1pe15.5))')
 891      &         vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
 892             enddo
 893           enddo
 894         enddo
 895       endif
 896 C
 897 C
 898 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
 899 C         interaction energy of the Gly, Ala, and Pro prototypes.
 900 C
 901       if (lprint) then
 902         write (iout,*)
 903         write (iout,*) "Coefficients of the cumulants"
 904       endif
 905       read (ifourier,*) nloctyp
 906       do i=0,nloctyp-1
 907         read (ifourier,*,end=115,err=115)
 908         read (ifourier,*,end=115,err=115) (b(ii),ii=1,13)
 909         if (lprint) then
 910         write (iout,*) 'Type',i
 911         write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
 912         endif
 913         B1(1,i)  = b(3)
 914         B1(2,i)  = b(5)
 915         B1(1,-i) = b(3)
 916         B1(2,-i) = -b(5)
 917 c        b1(1,i)=0.0d0
 918 c        b1(2,i)=0.0d0
 919         B1tilde(1,i) = b(3)
 920         B1tilde(2,i) =-b(5)
 921         B1tilde(1,-i) =-b(3)
 922         B1tilde(2,-i) =b(5)
 923 c        b1tilde(1,i)=0.0d0
 924 c        b1tilde(2,i)=0.0d0
 925         B2(1,i)  = b(2)
 926         B2(2,i)  = b(4)
 927         B2(1,-i)  =b(2)
 928         B2(2,-i)  =-b(4)
 929
 930 c        b2(1,i)=0.0d0
 931 c        b2(2,i)=0.0d0
 932         CC(1,1,i)= b(7)
 933         CC(2,2,i)=-b(7)
 934         CC(2,1,i)= b(9)
 935         CC(1,2,i)= b(9)
 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 c        CC(1,1,i)=0.0d0
 941 c        CC(2,2,i)=0.0d0
 942 c        CC(2,1,i)=0.0d0
 943 c        CC(1,2,i)=0.0d0
 944         Ctilde(1,1,i)=b(7)
 945         Ctilde(1,2,i)=b(9)
 946         Ctilde(2,1,i)=-b(9)
 947         Ctilde(2,2,i)=b(7)
 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
 953 c        Ctilde(1,1,i)=0.0d0
 954 c        Ctilde(1,2,i)=0.0d0
 955 c        Ctilde(2,1,i)=0.0d0
 956 c        Ctilde(2,2,i)=0.0d0
 957         DD(1,1,i)= b(6)
 958         DD(2,2,i)=-b(6)
 959         DD(2,1,i)= b(8)
 960         DD(1,2,i)= b(8)
 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 c        DD(1,1,i)=0.0d0
 966 c        DD(2,2,i)=0.0d0
 967 c        DD(2,1,i)=0.0d0
 968 c        DD(1,2,i)=0.0d0
 969         Dtilde(1,1,i)=b(6)
 970         Dtilde(1,2,i)=b(8)
 971         Dtilde(2,1,i)=-b(8)
 972         Dtilde(2,2,i)=b(6)
 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
 978 c        Dtilde(1,1,i)=0.0d0
 979 c        Dtilde(1,2,i)=0.0d0
 980 c        Dtilde(2,1,i)=0.0d0
 981 c        Dtilde(2,2,i)=0.0d0
 982         EE(1,1,i)= b(10)+b(11)
 983         EE(2,2,i)=-b(10)+b(11)
 984         EE(2,1,i)= b(12)-b(13)
 985         EE(1,2,i)= b(12)+b(13)
 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
 991 c        ee(1,1,i)=1.0d0
 992 c        ee(2,2,i)=1.0d0
 993 c        ee(2,1,i)=0.0d0
 994 c        ee(1,2,i)=0.0d0
 995 c        ee(2,1,i)=ee(1,2,i)
 996       enddo
 997       if (lprint) then
 998       do i=1,nloctyp
 999         write (iout,*) 'Type',i
1000         write (iout,*) 'B1'
1001         write(iout,*) B1(1,i),B1(2,i)
1002         write (iout,*) 'B2'
1003         write(iout,*) B2(1,i),B2(2,i)
1004         write (iout,*) 'CC'
1005         do j=1,2
1006           write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
1007         enddo
1008         write(iout,*) 'DD'
1009         do j=1,2
1010           write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
1011         enddo
1012         write(iout,*) 'EE'
1013         do j=1,2
1014           write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
1015         enddo
1016       enddo
1017       endif
1018 C
1019 C Read electrostatic-interaction parameters
1020 C
1021       if (lprint) then
1022         write (iout,*)
1023         write (iout,'(/a)') 'Electrostatic interaction constants:'
1024         write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
1025      &            'IT','JT','APP','BPP','AEL6','AEL3'
1026       endif
1027       read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
1028       read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
1029       read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
1030       read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
1031       close (ielep)
1032       do i=1,2
1033         do j=1,2
1034         rri=rpp(i,j)**6
1035         app (i,j)=epp(i,j)*rri*rri
1036         bpp (i,j)=-2.0D0*epp(i,j)*rri
1037         ael6(i,j)=elpp6(i,j)*4.2D0**6
1038         ael3(i,j)=elpp3(i,j)*4.2D0**3
1039         if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
1040      &                    ael6(i,j),ael3(i,j)
1041         enddo
1042       enddo
1043 C
1044 C Read side-chain interaction parameters.
1045 C
1046       read (isidep,*,end=117,err=117) ipot,expon
1047       if (ipot.lt.1 .or. ipot.gt.5) then
1048         write (iout,'(2a)') 'Error while reading SC interaction',
1049      &               'potential file - unknown potential type.'
1050 #ifdef MPI
1051         call MPI_Finalize(Ierror)
1052 #endif
1053         stop
1054       endif
1055       expon2=expon/2
1056       if(me.eq.king)
1057      & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
1058      & ', exponents are ',expon,2*expon
1059       goto (10,20,30,30,40) ipot
1060 C----------------------- LJ potential ---------------------------------
1061    10 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
1062      &   (sigma0(i),i=1,ntyp)
1063       if (lprint) then
1064         write (iout,'(/a/)') 'Parameters of the LJ potential:'
1065         write (iout,'(a/)') 'The epsilon array:'
1066         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1067         write (iout,'(/a)') 'One-body parameters:'
1068         write (iout,'(a,4x,a)') 'residue','sigma'
1069         write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1070       endif
1071       goto 50
1072 C----------------------- LJK potential --------------------------------
1073    20 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
1074      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1075       if (lprint) then
1076         write (iout,'(/a/)') 'Parameters of the LJK potential:'
1077         write (iout,'(a/)') 'The epsilon array:'
1078         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1079         write (iout,'(/a)') 'One-body parameters:'
1080         write (iout,'(a,4x,2a)') 'residue','   sigma  ','    r0    '
1081         write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
1082      &        i=1,ntyp)
1083       endif
1084       goto 50
1085 C---------------------- GB or BP potential -----------------------------
1086    30 do i=1,ntyp
1087         read (isidep,*,end=116,err=116)(eps(i,j),j=i,ntyp)
1088        enddo
1089        read (isidep,*,end=116,err=116)(sigma0(i),i=1,ntyp)
1090        read (isidep,*,end=116,err=116)(sigii(i),i=1,ntyp)
1091        read (isidep,*,end=116,err=116)(chip(i),i=1,ntyp)
1092        read (isidep,*,end=116,err=116)(alp(i),i=1,ntyp)
1093
1094 c   30 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
1095 c     &  (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip(i),i=1,ntyp),
1096 c     &  (alp(i),i=1,ntyp)
1097 C For the GB potential convert sigma'**2 into chi'
1098       if (ipot.eq.4) then
1099         do i=1,ntyp
1100           chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1101         enddo
1102       endif
1103       if (lprint) then
1104         write (iout,'(/a/)') 'Parameters of the BP potential:'
1105         write (iout,'(a/)') 'The epsilon array:'
1106         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1107         write (iout,'(/a)') 'One-body parameters:'
1108         write (iout,'(a,4x,4a)') 'residue','   sigma  ','s||/s_|_^2',
1109      &       '    chip  ','    alph  '
1110         write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1111      &                     chip(i),alp(i),i=1,ntyp)
1112       endif
1113       goto 50
1114 C--------------------- GBV potential -----------------------------------
1115    40 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
1116      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1117      &  (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1118       if (lprint) then
1119         write (iout,'(/a/)') 'Parameters of the GBV potential:'
1120         write (iout,'(a/)') 'The epsilon array:'
1121         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1122         write (iout,'(/a)') 'One-body parameters:'
1123         write (iout,'(a,4x,5a)') 'residue','   sigma  ','    r0    ',
1124      &      's||/s_|_^2','    chip  ','    alph  '
1125         write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1126      &           sigii(i),chip(i),alp(i),i=1,ntyp)
1127       endif
1128    50 continue
1129       close (isidep)
1130 C-----------------------------------------------------------------------
1131 C Calculate the "working" parameters of SC interactions.
1132       do i=2,ntyp
1133         do j=1,i-1
1134           eps(i,j)=eps(j,i)
1135         enddo
1136       enddo
1137       do i=1,ntyp
1138         do j=i,ntyp
1139           sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1140           sigma(j,i)=sigma(i,j)
1141           rs0(i,j)=dwa16*sigma(i,j)
1142           rs0(j,i)=rs0(i,j)
1143         enddo
1144       enddo
1145       if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1146      & 'Working parameters of the SC interactions:',
1147      & '     a    ','     b    ','   augm   ','  sigma ','   r0   ',
1148      & '  chi1   ','   chi2   '
1149       do i=1,ntyp
1150         do j=i,ntyp
1151           epsij=eps(i,j)
1152           if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1153             rrij=sigma(i,j)
1154           else
1155             rrij=rr0(i)+rr0(j)
1156           endif
1157           r0(i,j)=rrij
1158           r0(j,i)=rrij
1159           rrij=rrij**expon
1160           epsij=eps(i,j)
1161           sigeps=dsign(1.0D0,epsij)
1162           epsij=dabs(epsij)
1163           aa(i,j)=epsij*rrij*rrij
1164           bb(i,j)=-sigeps*epsij*rrij
1165           aa(j,i)=aa(i,j)
1166           bb(j,i)=bb(i,j)
1167           if (ipot.gt.2) then
1168             sigt1sq=sigma0(i)**2
1169             sigt2sq=sigma0(j)**2
1170             sigii1=sigii(i)
1171             sigii2=sigii(j)
1172             ratsig1=sigt2sq/sigt1sq
1173             ratsig2=1.0D0/ratsig1
1174             chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1175             if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1176             rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1177           else
1178             rsum_max=sigma(i,j)
1179           endif
1180 c         if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1181             sigmaii(i,j)=rsum_max
1182             sigmaii(j,i)=rsum_max
1183 c         else
1184 c           sigmaii(i,j)=r0(i,j)
1185 c           sigmaii(j,i)=r0(i,j)
1186 c         endif
1187 cd        write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1188           if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1189             r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1190             augm(i,j)=epsij*r_augm**(2*expon)
1191 c           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1192             augm(j,i)=augm(i,j)
1193           else
1194             augm(i,j)=0.0D0
1195             augm(j,i)=0.0D0
1196           endif
1197           if (lprint) then
1198             write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1199      &      restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1200      &      sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1201           endif
1202         enddo
1203       enddo
1204 #ifdef OLDSCP
1205 C
1206 C Define the SC-p interaction constants (hard-coded; old style)
1207 C
1208       do i=1,ntyp
1209 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1210 C helix formation)
1211 c       aad(i,1)=0.3D0*4.0D0**12
1212 C Following line for constants currently implemented
1213 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1214         aad(i,1)=1.5D0*4.0D0**12
1215 c       aad(i,1)=0.17D0*5.6D0**12
1216         aad(i,2)=aad(i,1)
1217 C "Soft" SC-p repulsion
1218         bad(i,1)=0.0D0
1219 C Following line for constants currently implemented
1220 c       aad(i,1)=0.3D0*4.0D0**6
1221 C "Hard" SC-p repulsion
1222         bad(i,1)=3.0D0*4.0D0**6
1223 c       bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1224         bad(i,2)=bad(i,1)
1225 c       aad(i,1)=0.0D0
1226 c       aad(i,2)=0.0D0
1227 c       bad(i,1)=1228.8D0
1228 c       bad(i,2)=1228.8D0
1229       enddo
1230 #else
1231 C
1232 C 8/9/01 Read the SC-p interaction constants from file
1233 C
1234       do i=1,ntyp
1235         read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1236       enddo
1237       do i=1,ntyp
1238         aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1239         aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1240         bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1241         bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1242       enddo
1243
1244       if (lprint) then
1245         write (iout,*) "Parameters of SC-p interactions:"
1246         do i=1,ntyp
1247           write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1248      &     eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1249         enddo
1250       endif
1251 #endif
1252 C
1253 C Define the constants of the disulfide bridge
1254 C
1255       ebr=-5.50D0
1256 c
1257 c Old arbitrary potential - commented out.
1258 c
1259 c      dbr= 4.20D0
1260 c      fbr= 3.30D0
1261 c
1262 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1263 c energy surface of diethyl disulfide.
1264 c A. Liwo and U. Kozlowska, 11/24/03
1265 c
1266       D0CM = 3.78d0
1267       AKCM = 15.1d0
1268       AKTH = 11.0d0
1269       AKCT = 12.0d0
1270       V1SS =-1.08d0
1271       V2SS = 7.61d0
1272       V3SS = 13.7d0
1273 c      akcm=0.0d0
1274 c      akth=0.0d0
1275 c      akct=0.0d0
1276 c      v1ss=0.0d0
1277 c      v2ss=0.0d0
1278 c      v3ss=0.0d0
1279
1280       if(me.eq.king) then
1281       write (iout,'(/a)') "Disulfide bridge parameters:"
1282       write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1283       write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1284       write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1285       write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1286      &  ' v3ss:',v3ss
1287       endif
1288       return
1289   111 write (iout,*) "Error reading bending energy parameters."
1290       goto 999
1291   112 write (iout,*) "Error reading rotamer energy parameters."
1292       goto 999
1293   113 write (iout,*) "Error reading torsional energy parameters."
1294       goto 999
1295   114 write (iout,*) "Error reading double torsional energy parameters."
1296       goto 999
1297   115 write (iout,*)
1298      &  "Error reading cumulant (multibody energy) parameters."
1299       goto 999
1300   116 write (iout,*) "Error reading electrostatic energy parameters."
1301       goto 999
1302   117 write (iout,*) "Error reading side chain interaction parameters."
1303       goto 999
1304   118 write (iout,*) "Error reading SCp interaction parameters."
1305       goto 999
1306   119 write (iout,*) "Error reading SCCOR parameters"
1307   999 continue
1308 #ifdef MPI
1309       call MPI_Finalize(Ierror)
1310 #endif
1311       stop
1312       return
1313       end
1314
1315
1316       subroutine getenv_loc(var, val)
1317       character(*) var, val
1318
1319 #ifdef WINIFL
1320       character(2000) line
1321       external ilen
1322
1323       open (196,file='env',status='old',readonly,shared)
1324       iread=0
1325 c      write(*,*)'looking for ',var
1326 10    read(196,*,err=11,end=11)line
1327       iread=index(line,var)
1328 c      write(*,*)iread,' ',var,' ',line
1329       if (iread.eq.0) go to 10
1330 c      write(*,*)'---> ',line
1331 11    continue
1332       if(iread.eq.0) then
1333 c       write(*,*)'CHUJ'
1334        val=''
1335       else
1336        iread=iread+ilen(var)+1
1337        read (line(iread:),*,err=12,end=12) val
1338 c       write(*,*)'OK: ',var,' = ',val
1339       endif
1340       close(196)
1341       return
1342 12    val=''
1343       close(196)
1344 #elif (defined CRAY)
1345       integer lennam,lenval,ierror
1346 c
1347 c        getenv using a POSIX call, useful on the T3D
1348 c        Sept 1996, comment out error check on advice of H. Pritchard
1349 c
1350       lennam = len(var)
1351       if(lennam.le.0) stop '--error calling getenv--'
1352       call pxfgetenv(var,lennam,val,lenval,ierror)
1353 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1354 #else
1355       call getenv(var,val)
1356 #endif
1357
1358       return
1359       end