source/unres/src_MD-M/parmread.F

   1       subroutine parmread
   2 C
   3 C Read the parameters of the probability distributions of the virtual-bond
   4 C valence angles and the side chains and energy parameters.
   5 C
   6 C Important! Energy-term weights ARE NOT read here; they are read from the
   7 C main input file instead, because NO defaults have yet been set for these
   8 C parameters.
   9 C
  10       implicit real*8 (a-h,o-z)
  11       include 'DIMENSIONS'
  12 #ifdef MPI
  13       include "mpif.h"
  14       integer IERROR
  15 #endif
  16       include 'COMMON.IOUNITS'
  17       include 'COMMON.CHAIN'
  18       include 'COMMON.INTERACT'
  19       include 'COMMON.GEO'
  20       include 'COMMON.LOCAL'
  21       include 'COMMON.TORSION'
  22       include 'COMMON.SCCOR'
  23       include 'COMMON.SCROT'
  24       include 'COMMON.FFIELD'
  25       include 'COMMON.NAMES'
  26       include 'COMMON.SBRIDGE'
  27       include 'COMMON.MD'
  28       include 'COMMON.SETUP'
  29       character*1 t1,t2,t3
  30       character*1 onelett(4) /"G","A","P","D"/
  31       character*1 toronelet(-2:2) /"p","a","G","A","P"/
  32       logical lprint,LaTeX
  33       dimension blower(3,3,maxlob)
  34       dimension b(13)
  35       character*3 lancuch,ucase
  36 C
  37 C For printing parameters after they are read set the following in the UNRES
  38 C C-shell script:
  39 C
  40 C setenv PRINT_PARM YES
  41 C
  42 C To print parameters in LaTeX format rather than as ASCII tables:
  43 C
  44 C setenv LATEX YES
  45 C
  46       call getenv_loc("PRINT_PARM",lancuch)
  47       lprint = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
  48       call getenv_loc("LATEX",lancuch)
  49       LaTeX = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
  50 C
  51       dwa16=2.0d0**(1.0d0/6.0d0)
  52       itypro=20
  53 C Assign virtual-bond length
  54       vbl=3.8D0
  55       vblinv=1.0D0/vbl
  56       vblinv2=vblinv*vblinv
  57 c
  58 c Read the virtual-bond parameters, masses, and moments of inertia
  59 c and Stokes' radii of the peptide group and side chains
  60 c
  61 #ifdef CRYST_BOND
  62       read (ibond,*) vbldp0,akp,mp,ip,pstok
  63       do i=1,ntyp
  64         nbondterm(i)=1
  65         read (ibond,*) vbldsc0(1,i),aksc(1,i),msc(i),isc(i),restok(i)
  66         dsc(i) = vbldsc0(1,i)
  67         if (i.eq.10) then
  68           dsc_inv(i)=0.0D0
  69         else
  70           dsc_inv(i)=1.0D0/dsc(i)
  71         endif
  72       enddo
  73 #else
  74       read (ibond,*) junk,vbldp0,akp,rjunk,mp,ip,pstok
  75       do i=1,ntyp
  76         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"
 392       do i=1,ntyp
 393         read (ithep,*,err=111,end=111) a0thet(i),(athet(j,i,1,1),j=1,2),
 394      &    (bthet(j,i,1,1),j=1,2)
 395         read (ithep,*,err=111,end=111) (polthet(j,i),j=0,3)
 396         read (ithep,*,err=111,end=111) (gthet(j,i),j=1,3)
 397         read (ithep,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
 398         sigc0(i)=sigc0(i)**2
 399       enddo
 400       do i=1,ntyp
 401       athet(1,i,1,-1)=athet(1,i,1,1)
 402       athet(2,i,1,-1)=athet(2,i,1,1)
 403       bthet(1,i,1,-1)=-bthet(1,i,1,1)
 404       bthet(2,i,1,-1)=-bthet(2,i,1,1)
 405       athet(1,i,-1,1)=-athet(1,i,1,1)
 406       athet(2,i,-1,1)=-athet(2,i,1,1)
 407       bthet(1,i,-1,1)=bthet(1,i,1,1)
 408       bthet(2,i,-1,1)=bthet(2,i,1,1)
 409       enddo
 410       do i=-ntyp,-1
 411       a0thet(i)=a0thet(-i)
 412       athet(1,i,-1,-1)=athet(1,-i,1,1)
 413       athet(2,i,-1,-1)=-athet(2,-i,1,1)
 414       bthet(1,i,-1,-1)=bthet(1,-i,1,1)
 415       bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
 416       athet(1,i,-1,1)=athet(1,-i,1,1)
 417       athet(2,i,-1,1)=-athet(2,-i,1,1)
 418       bthet(1,i,-1,1)=-bthet(1,-i,1,1)
 419       bthet(2,i,-1,1)=bthet(2,-i,1,1)
 420       athet(1,i,1,-1)=-athet(1,-i,1,1)
 421       athet(2,i,1,-1)=athet(2,-i,1,1)
 422       bthet(1,i,1,-1)=bthet(1,-i,1,1)
 423       bthet(2,i,1,-1)=-bthet(2,-i,1,1)
 424       theta0(i)=theta0(-i)
 425       sig0(i)=sig0(-i)
 426       sigc0(i)=sigc0(-i)
 427        do j=0,3
 428         polthet(j,i)=polthet(j,-i)
 429        enddo
 430        do j=1,3
 431          gthet(j,i)=gthet(j,-i)
 432        enddo
 433       enddo
 434       write (2,*) "End reading THETA_PDB"
 435       close (ithep_pdb)
 436 #endif
 437       close(ithep)
 438 #ifdef CRYST_SC
 439 C
 440 C Read the parameters of the probability distribution/energy expression
 441 C of the side chains.
 442 C
 443       do i=1,ntyp
 444         read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
 445         if (i.eq.10) then
 446           dsc_inv(i)=0.0D0
 447         else
 448           dsc_inv(i)=1.0D0/dsc(i)
 449         endif
 450         if (i.ne.10) then
 451         do j=1,nlob(i)
 452           do k=1,3
 453             do l=1,3
 454               blower(l,k,j)=0.0D0
 455             enddo
 456           enddo
 457         enddo
 458         bsc(1,i)=0.0D0
 459         read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),
 460      &    ((blower(k,l,1),l=1,k),k=1,3)
 461         censc(1,1,-i)=censc(1,1,i)
 462         censc(2,1,-i)=censc(2,1,i)
 463         censc(3,1,-i)=-censc(3,1,i)
 464         do j=2,nlob(i)
 465           read (irotam,*,end=112,err=112) bsc(j,i)
 466           read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),
 467      &                                 ((blower(k,l,j),l=1,k),k=1,3)
 468         censc(1,j,-i)=censc(1,j,i)
 469         censc(2,j,-i)=censc(2,j,i)
 470         censc(3,j,-i)=-censc(3,j,i)
 471 C BSC is amplitude of Gaussian
 472         enddo
 473         do j=1,nlob(i)
 474           do k=1,3
 475             do l=1,k
 476               akl=0.0D0
 477               do m=1,3
 478                 akl=akl+blower(k,m,j)*blower(l,m,j)
 479               enddo
 480               gaussc(k,l,j,i)=akl
 481               gaussc(l,k,j,i)=akl
 482              if (((k.eq.3).and.(l.ne.3))
 483      &        .or.((l.eq.3).and.(k.ne.3))) then
 484                 gaussc(k,l,j,-i)=-akl
 485                 gaussc(l,k,j,-i)=-akl
 486               else
 487                 gaussc(k,l,j,-i)=akl
 488                 gaussc(l,k,j,-i)=akl
 489               endif
 490             enddo
 491           enddo
 492         enddo
 493         endif
 494       enddo
 495       close (irotam)
 496       if (lprint) then
 497         write (iout,'(/a)') 'Parameters of side-chain local geometry'
 498         do i=1,ntyp
 499           nlobi=nlob(i)
 500           if (nlobi.gt.0) then
 501             if (LaTeX) then
 502               write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
 503      &         ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
 504                write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
 505      &                             'log h',(bsc(j,i),j=1,nlobi)
 506                write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
 507      &        'x',((censc(k,j,i),k=1,3),j=1,nlobi)
 508               do k=1,3
 509                 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
 510      &                 ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
 511               enddo
 512             else
 513               write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
 514               write (iout,'(a,f10.4,4(16x,f10.4))')
 515      &                             'Center  ',(bsc(j,i),j=1,nlobi)
 516               write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),
 517      &           j=1,nlobi)
 518               write (iout,'(a)')
 519             endif
 520           endif
 521         enddo
 522       endif
 523 #else
 524 C
 525 C Read scrot parameters for potentials determined from all-atom AM1 calculations
 526 C added by Urszula Kozlowska 07/11/2007
 527 C
 528       do i=1,ntyp
 529         read (irotam,*,end=112,err=112)
 530        if (i.eq.10) then
 531          read (irotam,*,end=112,err=112)
 532        else
 533          do j=1,65
 534            read(irotam,*,end=112,err=112) sc_parmin(j,i)
 535          enddo
 536        endif
 537       enddo
 538 C
 539 C Read the parameters of the probability distribution/energy expression
 540 C of the side chains.
 541 C
 542       do i=1,ntyp
 543         read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
 544         if (i.eq.10) then
 545           dsc_inv(i)=0.0D0
 546         else
 547           dsc_inv(i)=1.0D0/dsc(i)
 548         endif
 549         if (i.ne.10) then
 550         do j=1,nlob(i)
 551           do k=1,3
 552             do l=1,3
 553               blower(l,k,j)=0.0D0
 554             enddo
 555           enddo
 556         enddo
 557         bsc(1,i)=0.0D0
 558         read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),
 559      &    ((blower(k,l,1),l=1,k),k=1,3)
 560         do j=2,nlob(i)
 561           read (irotam_pdb,*,end=112,err=112) bsc(j,i)
 562           read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),
 563      &                                 ((blower(k,l,j),l=1,k),k=1,3)
 564         enddo
 565         do j=1,nlob(i)
 566           do k=1,3
 567             do l=1,k
 568               akl=0.0D0
 569               do m=1,3
 570                 akl=akl+blower(k,m,j)*blower(l,m,j)
 571               enddo
 572               gaussc(k,l,j,i)=akl
 573               gaussc(l,k,j,i)=akl
 574             enddo
 575           enddo
 576         enddo
 577         endif
 578       enddo
 579       close (irotam_pdb)
 580 #endif
 581       close(irotam)
 582
 583 #ifdef CRYST_TOR
 584 C
 585 C Read torsional parameters in old format
 586 C
 587       read (itorp,*,end=113,err=113) ntortyp,nterm_old
 588       if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
 589       read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
 590       do i=1,ntortyp
 591         do j=1,ntortyp
 592           read (itorp,'(a)')
 593           do k=1,nterm_old
 594             read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
 595           enddo
 596         enddo
 597       enddo
 598       close (itorp)
 599       if (lprint) then
 600         write (iout,'(/a/)') 'Torsional constants:'
 601         do i=1,ntortyp
 602           do j=1,ntortyp
 603             write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
 604             write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
 605           enddo
 606         enddo
 607       endif
 608 #else
 609 C
 610 C Read torsional parameters
 611 C
 612       read (itorp,*,end=113,err=113) ntortyp
 613       read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
 614       do iblock=1,2
 615       do i=-ntyp,-1
 616        itortyp(i)=-itortyp(-i)
 617       enddo
 618       write (iout,*) 'ntortyp',ntortyp
 619       do i=0,ntortyp-1
 620         do j=-ntortyp+1,ntortyp-1
 621           read (itorp,*,end=113,err=113) nterm(i,j,iblock),
 622      &          nlor(i,j,iblock)
 623           nterm(-i,-j,iblock)=nterm(i,j,iblock)
 624           nlor(-i,-j,iblock)=nlor(i,j,iblock)
 625           v0ij=0.0d0
 626           si=-1.0d0
 627           do k=1,nterm(i,j,iblock)
 628             read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),
 629      &      v2(k,i,j,iblock)
 630             v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
 631             v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
 632             v0ij=v0ij+si*v1(k,i,j,iblock)
 633             si=-si
 634 c         write(iout,*) i,j,k,iblock,nterm(i,j,iblock)
 635 c         write(iout,*) v1(k,-i,-j,iblock),v1(k,i,j,iblock),
 636 c      &v2(k,-i,-j,iblock),v2(k,i,j,iblock)
 637           enddo
 638           do k=1,nlor(i,j,iblock)
 639             read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),
 640      &        vlor2(k,i,j),vlor3(k,i,j)
 641             v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
 642           enddo
 643           v0(i,j,iblock)=v0ij
 644           v0(-i,-j,iblock)=v0ij
 645         enddo
 646       enddo
 647       enddo
 648       close (itorp)
 649       if (lprint) then
 650         write (iout,'(/a/)') 'Torsional constants:'
 651         do i=1,ntortyp
 652           do j=1,ntortyp
 653             write (iout,*) 'ityp',i,' jtyp',j
 654             write (iout,*) 'Fourier constants'
 655             do k=1,nterm(i,j,iblock)
 656               write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
 657      &        v2(k,i,j,iblock)
 658             enddo
 659             write (iout,*) 'Lorenz constants'
 660             do k=1,nlor(i,j,iblock)
 661               write (iout,'(3(1pe15.5))')
 662      &         vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
 663             enddo
 664           enddo
 665         enddo
 666       endif
 667
 668 C
 669 C 6/23/01 Read parameters for double torsionals
 670 C
 671       do iblock=1,2
 672       do i=0,ntortyp-1
 673         do j=-ntortyp+1,ntortyp-1
 674           do k=-ntortyp+1,ntortyp-1
 675             read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
 676 c              write (iout,*) "OK onelett",
 677 c     &         i,j,k,t1,t2,t3
 678
 679             if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
 680      &        .or. t3.ne.toronelet(k)) then
 681               write (iout,*) "Error in double torsional parameter file",
 682      &         i,j,k,t1,t2,t3
 683 #ifdef MPI
 684               call MPI_Finalize(Ierror)
 685 #endif
 686                stop "Error in double torsional parameter file"
 687             endif
 688            read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),
 689      &         ntermd_2(i,j,k,iblock)
 690             ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
 691             ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
 692             read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,
 693      &         ntermd_1(i,j,k,iblock))
 694             read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,
 695      &         ntermd_1(i,j,k,iblock))
 696             read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,
 697      &         ntermd_1(i,j,k,iblock))
 698             read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,
 699      &         ntermd_1(i,j,k,iblock))
 700 C Martix of D parameters for one dimesional foureir series
 701             do l=1,ntermd_1(i,j,k,iblock)
 702              v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
 703              v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
 704              v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
 705              v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
 706 c            write(iout,*) "whcodze" ,
 707 c     & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
 708             enddo
 709             read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),
 710      &         v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
 711      &         v2s(m,l,i,j,k,iblock),
 712      &         m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
 713 C Martix of D parameters for two dimesional fourier series
 714             do l=1,ntermd_2(i,j,k,iblock)
 715              do m=1,l-1
 716              v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
 717              v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
 718              v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
 719              v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
 720              enddo!m
 721             enddo!l
 722           enddo!k
 723         enddo!j
 724       enddo!i
 725       enddo!iblock
 726       if (lprint) then
 727       write (iout,*)
 728       write (iout,*) 'Constants for double torsionals'
 729       do iblock=1,2
 730       do i=0,ntortyp-1
 731         do j=-ntortyp+1,ntortyp-1
 732           do k=-ntortyp+1,ntortyp-1
 733             write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
 734      &        ' nsingle',ntermd_1(i,j,k,iblock),
 735      &        ' ndouble',ntermd_2(i,j,k,iblock)
 736             write (iout,*)
 737             write (iout,*) 'Single angles:'
 738             do l=1,ntermd_1(i,j,k,iblock)
 739               write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
 740      &           v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
 741      &           v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),
 742      &           v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
 743             enddo
 744             write (iout,*)
 745             write (iout,*) 'Pairs of angles:'
 746             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
 747             do l=1,ntermd_2(i,j,k,iblock)
 748               write (iout,'(i5,20f10.5)')
 749      &         l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
 750             enddo
 751             write (iout,*)
 752             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
 753             do l=1,ntermd_2(i,j,k,iblock)
 754               write (iout,'(i5,20f10.5)')
 755      &         l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),
 756      &         (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
 757             enddo
 758             write (iout,*)
 759           enddo
 760         enddo
 761       enddo
 762       enddo
 763       endif
 764 #endif
 765 C Read of Side-chain backbone correlation parameters
 766 C Modified 11 May 2012 by Adasko
 767 CCC
 768 C
 769       read (isccor,*,end=119,err=119) nsccortyp
 770 #ifdef SCCORPDB
 771       read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
 772       do i=-ntyp,-1
 773         isccortyp(i)=-isccortyp(-i)
 774       enddo
 775       iscprol=isccortyp(20)
 776 c      write (iout,*) 'ntortyp',ntortyp
 777       maxinter=3
 778 cc maxinter is maximum interaction sites
 779       do l=1,maxinter
 780       do i=1,nsccortyp
 781         do j=1,nsccortyp
 782           read (isccor,*,end=119,err=119)
 783      &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
 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
1019 C
1020 C Read electrostatic-interaction parameters
1021 C
1022       if (lprint) then
1023         write (iout,*)
1024         write (iout,'(/a)') 'Electrostatic interaction constants:'
1025         write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
1026      &            'IT','JT','APP','BPP','AEL6','AEL3'
1027       endif
1028       read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
1029       read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
1030       read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
1031       read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
1032       close (ielep)
1033       do i=1,2
1034         do j=1,2
1035         rri=rpp(i,j)**6
1036         app (i,j)=epp(i,j)*rri*rri
1037         bpp (i,j)=-2.0D0*epp(i,j)*rri
1038         ael6(i,j)=elpp6(i,j)*4.2D0**6
1039         ael3(i,j)=elpp3(i,j)*4.2D0**3
1040 c        lprint=.true.
1041         if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
1042      &                    ael6(i,j),ael3(i,j)
1043 c        lprint=.false.
1044         enddo
1045       enddo
1046 C
1047 C Read side-chain interaction parameters.
1048 C
1049       read (isidep,*,end=117,err=117) ipot,expon
1050       if (ipot.lt.1 .or. ipot.gt.5) then
1051         write (iout,'(2a)') 'Error while reading SC interaction',
1052      &               'potential file - unknown potential type.'
1053 #ifdef MPI
1054         call MPI_Finalize(Ierror)
1055 #endif
1056         stop
1057       endif
1058       expon2=expon/2
1059       if(me.eq.king)
1060      & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
1061      & ', exponents are ',expon,2*expon
1062       goto (10,20,30,30,40) ipot
1063 C----------------------- LJ potential ---------------------------------
1064    10 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
1065      &   (sigma0(i),i=1,ntyp)
1066       if (lprint) then
1067         write (iout,'(/a/)') 'Parameters of the LJ potential:'
1068         write (iout,'(a/)') 'The epsilon array:'
1069         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1070         write (iout,'(/a)') 'One-body parameters:'
1071         write (iout,'(a,4x,a)') 'residue','sigma'
1072         write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1073       endif
1074       goto 50
1075 C----------------------- LJK potential --------------------------------
1076    20 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
1077      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1078       if (lprint) then
1079         write (iout,'(/a/)') 'Parameters of the LJK potential:'
1080         write (iout,'(a/)') 'The epsilon array:'
1081         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1082         write (iout,'(/a)') 'One-body parameters:'
1083         write (iout,'(a,4x,2a)') 'residue','   sigma  ','    r0    '
1084         write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
1085      &        i=1,ntyp)
1086       endif
1087       goto 50
1088 C---------------------- GB or BP potential -----------------------------
1089    30 do i=1,ntyp
1090        read (isidep,*,end=116,err=116)(eps(i,j),j=i,ntyp)
1091       enddo
1092       read (isidep,*,end=116,err=116)(sigma0(i),i=1,ntyp)
1093       read (isidep,*,end=116,err=116)(sigii(i),i=1,ntyp)
1094       read (isidep,*,end=116,err=116)(chip(i),i=1,ntyp)
1095       read (isidep,*,end=116,err=116)(alp(i),i=1,ntyp)
1096 C For the GB potential convert sigma'**2 into chi'
1097       if (ipot.eq.4) then
1098         do i=1,ntyp
1099           chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1100         enddo
1101       endif
1102       if (lprint) then
1103         write (iout,'(/a/)') 'Parameters of the BP potential:'
1104         write (iout,'(a/)') 'The epsilon array:'
1105         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1106         write (iout,'(/a)') 'One-body parameters:'
1107         write (iout,'(a,4x,4a)') 'residue','   sigma  ','s||/s_|_^2',
1108      &       '    chip  ','    alph  '
1109         write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1110      &                     chip(i),alp(i),i=1,ntyp)
1111       endif
1112       goto 50
1113 C--------------------- GBV potential -----------------------------------
1114    40 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
1115      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1116      &  (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1117       if (lprint) then
1118         write (iout,'(/a/)') 'Parameters of the GBV potential:'
1119         write (iout,'(a/)') 'The epsilon array:'
1120         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1121         write (iout,'(/a)') 'One-body parameters:'
1122         write (iout,'(a,4x,5a)') 'residue','   sigma  ','    r0    ',
1123      &      's||/s_|_^2','    chip  ','    alph  '
1124         write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1125      &           sigii(i),chip(i),alp(i),i=1,ntyp)
1126       endif
1127    50 continue
1128       close (isidep)
1129 C-----------------------------------------------------------------------
1130 C Calculate the "working" parameters of SC interactions.
1131       do i=2,ntyp
1132         do j=1,i-1
1133           eps(i,j)=eps(j,i)
1134         enddo
1135       enddo
1136       do i=1,ntyp
1137         do j=i,ntyp
1138           sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1139           sigma(j,i)=sigma(i,j)
1140           rs0(i,j)=dwa16*sigma(i,j)
1141           rs0(j,i)=rs0(i,j)
1142         enddo
1143       enddo
1144       if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1145      & 'Working parameters of the SC interactions:',
1146      & '     a    ','     b    ','   augm   ','  sigma ','   r0   ',
1147      & '  chi1   ','   chi2   '
1148       do i=1,ntyp
1149         do j=i,ntyp
1150           epsij=eps(i,j)
1151           if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1152             rrij=sigma(i,j)
1153           else
1154             rrij=rr0(i)+rr0(j)
1155           endif
1156           r0(i,j)=rrij
1157           r0(j,i)=rrij
1158           rrij=rrij**expon
1159           epsij=eps(i,j)
1160           sigeps=dsign(1.0D0,epsij)
1161           epsij=dabs(epsij)
1162           aa(i,j)=epsij*rrij*rrij
1163           bb(i,j)=-sigeps*epsij*rrij
1164           aa(j,i)=aa(i,j)
1165           bb(j,i)=bb(i,j)
1166           if (ipot.gt.2) then
1167             sigt1sq=sigma0(i)**2
1168             sigt2sq=sigma0(j)**2
1169             sigii1=sigii(i)
1170             sigii2=sigii(j)
1171             ratsig1=sigt2sq/sigt1sq
1172             ratsig2=1.0D0/ratsig1
1173             chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1174             if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1175             rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1176           else
1177             rsum_max=sigma(i,j)
1178           endif
1179 c         if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1180             sigmaii(i,j)=rsum_max
1181             sigmaii(j,i)=rsum_max
1182 c         else
1183 c           sigmaii(i,j)=r0(i,j)
1184 c           sigmaii(j,i)=r0(i,j)
1185 c         endif
1186 cd        write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1187           if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1188             r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1189             augm(i,j)=epsij*r_augm**(2*expon)
1190 c           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1191             augm(j,i)=augm(i,j)
1192           else
1193             augm(i,j)=0.0D0
1194             augm(j,i)=0.0D0
1195           endif
1196           if (lprint) then
1197             write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1198      &      restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1199      &      sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1200           endif
1201         enddo
1202       enddo
1203 #ifdef OLDSCP
1204 C
1205 C Define the SC-p interaction constants (hard-coded; old style)
1206 C
1207       do i=1,ntyp
1208 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1209 C helix formation)
1210 c       aad(i,1)=0.3D0*4.0D0**12
1211 C Following line for constants currently implemented
1212 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1213         aad(i,1)=1.5D0*4.0D0**12
1214 c       aad(i,1)=0.17D0*5.6D0**12
1215         aad(i,2)=aad(i,1)
1216 C "Soft" SC-p repulsion
1217         bad(i,1)=0.0D0
1218 C Following line for constants currently implemented
1219 c       aad(i,1)=0.3D0*4.0D0**6
1220 C "Hard" SC-p repulsion
1221         bad(i,1)=3.0D0*4.0D0**6
1222 c       bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1223         bad(i,2)=bad(i,1)
1224 c       aad(i,1)=0.0D0
1225 c       aad(i,2)=0.0D0
1226 c       bad(i,1)=1228.8D0
1227 c       bad(i,2)=1228.8D0
1228       enddo
1229 #else
1230 C
1231 C 8/9/01 Read the SC-p interaction constants from file
1232 C
1233       do i=1,ntyp
1234         read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1235       enddo
1236       do i=1,ntyp
1237         aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1238         aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1239         bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1240         bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1241       enddo
1242 c      lprint=.true.
1243       if (lprint) then
1244         write (iout,*) "Parameters of SC-p interactions:"
1245         do i=1,ntyp
1246           write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1247      &     eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1248         enddo
1249       endif
1250 c      lprint=.false.
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