source/unres/src_MD/parmread.F

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