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