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       logical lprint,LaTeX
  32       dimension blower(3,3,maxlob)
  33       dimension b(13)
  34       character*3 lancuch,ucase
  35 C
  36 C For printing parameters after they are read set the following in the UNRES
  37 C C-shell script:
  38 C
  39 C setenv PRINT_PARM YES
  40 C
  41 C To print parameters in LaTeX format rather than as ASCII tables:
  42 C
  43 C setenv LATEX YES
  44 C
  45       call getenv_loc("PRINT_PARM",lancuch)
  46       lprint = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
  47       call getenv_loc("LATEX",lancuch)
  48       LaTeX = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
  49 C
  50       dwa16=2.0d0**(1.0d0/6.0d0)
  51       itypro=20
  52 C Assign virtual-bond length
  53       vbl=3.8D0
  54       vblinv=1.0D0/vbl
  55       vblinv2=vblinv*vblinv
  56 c
  57 c Read the virtual-bond parameters, masses, and moments of inertia
  58 c and Stokes' radii of the peptide group and side chains
  59 c
  60 #ifdef CRYST_BOND
  61       read (ibond,*) vbldp0,akp,mp,ip,pstok
  62       do i=1,ntyp
  63         nbondterm(i)=1
  64         read (ibond,*) vbldsc0(1,i),aksc(1,i),msc(i),isc(i),restok(i)
  65         dsc(i) = vbldsc0(1,i)
  66         if (i.eq.10) then
  67           dsc_inv(i)=0.0D0
  68         else
  69           dsc_inv(i)=1.0D0/dsc(i)
  70         endif
  71       enddo
  72 #else
  73       read (ibond,*) junk,vbldp0,akp,rjunk,mp,ip,pstok
  74       do i=1,ntyp
  75         read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),
  76      &   j=1,nbondterm(i)),msc(i),isc(i),restok(i)
  77         dsc(i) = vbldsc0(1,i)
  78         if (i.eq.10) then
  79           dsc_inv(i)=0.0D0
  80         else
  81           dsc_inv(i)=1.0D0/dsc(i)
  82         endif
  83       enddo
  84 #endif
  85       if (lprint) then
  86         write(iout,'(/a/)')"Dynamic constants of the interaction sites:"
  87         write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K','A0','mass',
  88      &   'inertia','Pstok'
  89         write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0,mp,ip,pstok
  90         do i=1,ntyp
  91           write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),
  92      &      vbldsc0(1,i),aksc(1,i),abond0(1,i),msc(i),isc(i),restok(i)
  93           do j=2,nbondterm(i)
  94             write (iout,'(13x,3f10.5)')
  95      &        vbldsc0(j,i),aksc(j,i),abond0(j,i)
  96           enddo
  97         enddo
  98       endif
  99 #ifdef CRYST_THETA
 100 C
 101 C Read the parameters of the probability distribution/energy expression
 102 C of the virtual-bond valence angles theta
 103 C
 104       do i=1,ntyp
 105         read (ithep,*,err=111,end=111) a0thet(i),(athet(j,i),j=1,2),
 106      &    (bthet(j,i),j=1,2)
 107         read (ithep,*,err=111,end=111) (polthet(j,i),j=0,3)
 108         read (ithep,*,err=111,end=111) (gthet(j,i),j=1,3)
 109         read (ithep,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
 110         sigc0(i)=sigc0(i)**2
 111       enddo
 112       close (ithep)
 113       if (lprint) then
 114       if (.not.LaTeX) then
 115         write (iout,'(a)')
 116      &    'Parameters of the virtual-bond valence angles:'
 117         write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
 118      & '    ATHETA0   ','         A1   ','        A2    ',
 119      & '        B1    ','         B2   '
 120         do i=1,ntyp
 121           write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
 122      &        a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2)
 123         enddo
 124         write (iout,'(/a/9x,5a/79(1h-))')
 125      & 'Parameters of the expression for sigma(theta_c):',
 126      & '     ALPH0    ','      ALPH1   ','     ALPH2    ',
 127      & '     ALPH3    ','    SIGMA0C   '
 128         do i=1,ntyp
 129           write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
 130      &      (polthet(j,i),j=0,3),sigc0(i)
 131         enddo
 132         write (iout,'(/a/9x,5a/79(1h-))')
 133      & 'Parameters of the second gaussian:',
 134      & '    THETA0    ','     SIGMA0   ','        G1    ',
 135      & '        G2    ','         G3   '
 136         do i=1,ntyp
 137           write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
 138      &       sig0(i),(gthet(j,i),j=1,3)
 139         enddo
 140        else
 141         write (iout,'(a)')
 142      &    'Parameters of the virtual-bond valence angles:'
 143         write (iout,'(/a/9x,5a/79(1h-))')
 144      & 'Coefficients of expansion',
 145      & '     theta0   ','    a1*10^2   ','   a2*10^2    ',
 146      & '   b1*10^1    ','    b2*10^1   '
 147         do i=1,ntyp
 148           write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),
 149      &        a0thet(i),(100*athet(j,i),j=1,2),(10*bthet(j,i),j=1,2)
 150         enddo
 151         write (iout,'(/a/9x,5a/79(1h-))')
 152      & 'Parameters of the expression for sigma(theta_c):',
 153      & ' alpha0       ','  alph1       ',' alph2        ',
 154      & ' alhp3        ','   sigma0c    '
 155         do i=1,ntyp
 156           write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),
 157      &      (polthet(j,i),j=0,3),sigc0(i)
 158         enddo
 159         write (iout,'(/a/9x,5a/79(1h-))')
 160      & 'Parameters of the second gaussian:',
 161      & '    theta0    ','  sigma0*10^2 ','      G1*10^-1',
 162      & '        G2    ','   G3*10^1    '
 163         do i=1,ntyp
 164           write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),
 165      &       100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
 166         enddo
 167       endif
 168       endif
 169 #else
 170 C
 171 C Read the parameters of Utheta determined from ab initio surfaces
 172 C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
 173 C
 174       read (ithep,*,err=111,end=111) nthetyp,ntheterm,ntheterm2,
 175      &  ntheterm3,nsingle,ndouble
 176       nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
 177       read (ithep,*,err=111,end=111) (ithetyp(i),i=1,ntyp1)
 178       do i=1,maxthetyp
 179         do j=1,maxthetyp
 180           do k=1,maxthetyp
 181             aa0thet(i,j,k)=0.0d0
 182             do l=1,ntheterm
 183               aathet(l,i,j,k)=0.0d0
 184             enddo
 185             do l=1,ntheterm2
 186               do m=1,nsingle
 187                 bbthet(m,l,i,j,k)=0.0d0
 188                 ccthet(m,l,i,j,k)=0.0d0
 189                 ddthet(m,l,i,j,k)=0.0d0
 190                 eethet(m,l,i,j,k)=0.0d0
 191               enddo
 192             enddo
 193             do l=1,ntheterm3
 194               do m=1,ndouble
 195                 do mm=1,ndouble
 196                  ffthet(mm,m,l,i,j,k)=0.0d0
 197                  ggthet(mm,m,l,i,j,k)=0.0d0
 198                 enddo
 199               enddo
 200             enddo
 201           enddo
 202         enddo
 203       enddo
 204       do i=1,nthetyp
 205         do j=1,nthetyp
 206           do k=1,nthetyp
 207             read (ithep,'(3a)',end=111,err=111) res1,res2,res3
 208             read (ithep,*,end=111,err=111) aa0thet(i,j,k)
 209             read (ithep,*,end=111,err=111)(aathet(l,i,j,k),l=1,ntheterm)
 210             read (ithep,*,end=111,err=111)
 211      &       ((bbthet(lll,ll,i,j,k),lll=1,nsingle),
 212      &        (ccthet(lll,ll,i,j,k),lll=1,nsingle),
 213      &        (ddthet(lll,ll,i,j,k),lll=1,nsingle),
 214      &        (eethet(lll,ll,i,j,k),lll=1,nsingle),ll=1,ntheterm2)
 215             read (ithep,*,end=111,err=111)
 216      &      (((ffthet(llll,lll,ll,i,j,k),ffthet(lll,llll,ll,i,j,k),
 217      &         ggthet(llll,lll,ll,i,j,k),ggthet(lll,llll,ll,i,j,k),
 218      &         llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
 219           enddo
 220         enddo
 221       enddo
 222 C
 223 C For dummy ends assign glycine-type coefficients of theta-only terms; the
 224 C coefficients of theta-and-gamma-dependent terms are zero.
 225 C
 226       do i=1,nthetyp
 227         do j=1,nthetyp
 228           do l=1,ntheterm
 229             aathet(l,i,j,nthetyp+1)=aathet(l,i,j,1)
 230             aathet(l,nthetyp+1,i,j)=aathet(l,1,i,j)
 231           enddo
 232           aa0thet(i,j,nthetyp+1)=aa0thet(i,j,1)
 233           aa0thet(nthetyp+1,i,j)=aa0thet(1,i,j)
 234         enddo
 235         do l=1,ntheterm
 236           aathet(l,nthetyp+1,i,nthetyp+1)=aathet(l,1,i,1)
 237         enddo
 238         aa0thet(nthetyp+1,i,nthetyp+1)=aa0thet(1,i,1)
 239       enddo
 240 C
 241 C Control printout of the coefficients of virtual-bond-angle potentials
 242 C
 243       if (lprint) then
 244         write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
 245         do i=1,nthetyp+1
 246           do j=1,nthetyp+1
 247             do k=1,nthetyp+1
 248               write (iout,'(//4a)')
 249      &         'Type ',onelett(i),onelett(j),onelett(k)
 250               write (iout,'(//a,10x,a)') " l","a[l]"
 251               write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k)
 252               write (iout,'(i2,1pe15.5)')
 253      &           (l,aathet(l,i,j,k),l=1,ntheterm)
 254             do l=1,ntheterm2
 255               write (iout,'(//2h m,4(9x,a,3h[m,,i1,1h]))')
 256      &          "b",l,"c",l,"d",l,"e",l
 257               do m=1,nsingle
 258                 write (iout,'(i2,4(1pe15.5))') m,
 259      &          bbthet(m,l,i,j,k),ccthet(m,l,i,j,k),
 260      &          ddthet(m,l,i,j,k),eethet(m,l,i,j,k)
 261               enddo
 262             enddo
 263             do l=1,ntheterm3
 264               write (iout,'(//3hm,n,4(6x,a,5h[m,n,,i1,1h]))')
 265      &          "f+",l,"f-",l,"g+",l,"g-",l
 266               do m=2,ndouble
 267                 do n=1,m-1
 268                   write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
 269      &              ffthet(n,m,l,i,j,k),ffthet(m,n,l,i,j,k),
 270      &              ggthet(n,m,l,i,j,k),ggthet(m,n,l,i,j,k)
 271                 enddo
 272               enddo
 273             enddo
 274           enddo
 275         enddo
 276       enddo
 277       call flush(iout)
 278       endif
 279 #endif
 280       close(ithep)
 281 #ifdef CRYST_SC
 282 C
 283 C Read the parameters of the probability distribution/energy expression
 284 C of the side chains.
 285 C
 286       do i=1,ntyp
 287         read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
 288         if (i.eq.10) then
 289           dsc_inv(i)=0.0D0
 290         else
 291           dsc_inv(i)=1.0D0/dsc(i)
 292         endif
 293         if (i.ne.10) then
 294         do j=1,nlob(i)
 295           do k=1,3
 296             do l=1,3
 297               blower(l,k,j)=0.0D0
 298             enddo
 299           enddo
 300         enddo
 301         bsc(1,i)=0.0D0
 302         read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),
 303      &    ((blower(k,l,1),l=1,k),k=1,3)
 304         do j=2,nlob(i)
 305           read (irotam,*,end=112,err=112) bsc(j,i)
 306           read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),
 307      &                                 ((blower(k,l,j),l=1,k),k=1,3)
 308         enddo
 309         do j=1,nlob(i)
 310           do k=1,3
 311             do l=1,k
 312               akl=0.0D0
 313               do m=1,3
 314                 akl=akl+blower(k,m,j)*blower(l,m,j)
 315               enddo
 316               gaussc(k,l,j,i)=akl
 317               gaussc(l,k,j,i)=akl
 318             enddo
 319           enddo
 320         enddo
 321         endif
 322       enddo
 323       close (irotam)
 324       if (lprint) then
 325         write (iout,'(/a)') 'Parameters of side-chain local geometry'
 326         do i=1,ntyp
 327           nlobi=nlob(i)
 328           if (nlobi.gt.0) then
 329             if (LaTeX) then
 330               write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
 331      &         ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
 332                write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
 333      &                             'log h',(bsc(j,i),j=1,nlobi)
 334                write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
 335      &        'x',((censc(k,j,i),k=1,3),j=1,nlobi)
 336               do k=1,3
 337                 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
 338      &                 ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
 339               enddo
 340             else
 341               write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
 342               write (iout,'(a,f10.4,4(16x,f10.4))')
 343      &                             'Center  ',(bsc(j,i),j=1,nlobi)
 344               write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),
 345      &           j=1,nlobi)
 346               write (iout,'(a)')
 347             endif
 348           endif
 349         enddo
 350       endif
 351 #else
 352 C
 353 C Read scrot parameters for potentials determined from all-atom AM1 calculations
 354 C added by Urszula Kozlowska 07/11/2007
 355 C
 356       do i=1,ntyp
 357         read (irotam,*,end=112,err=112)
 358        if (i.eq.10) then
 359          read (irotam,*,end=112,err=112)
 360        else
 361          do j=1,65
 362            read(irotam,*,end=112,err=112) sc_parmin(j,i)
 363          enddo
 364        endif
 365       enddo
 366 #endif
 367       close(irotam)
 368
 369 #ifdef CRYST_TOR
 370 C
 371 C Read torsional parameters in old format
 372 C
 373       read (itorp,*,end=113,err=113) ntortyp,nterm_old
 374       if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
 375       read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
 376       do i=1,ntortyp
 377         do j=1,ntortyp
 378           read (itorp,'(a)')
 379           do k=1,nterm_old
 380             read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
 381           enddo
 382         enddo
 383       enddo
 384       close (itorp)
 385       if (lprint) then
 386         write (iout,'(/a/)') 'Torsional constants:'
 387         do i=1,ntortyp
 388           do j=1,ntortyp
 389             write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
 390             write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
 391           enddo
 392         enddo
 393       endif
 394 #else
 395 C
 396 C Read torsional parameters
 397 C
 398       read (itorp,*,end=113,err=113) ntortyp
 399       read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
 400 c      write (iout,*) 'ntortyp',ntortyp
 401       do i=1,ntortyp
 402         do j=1,ntortyp
 403           read (itorp,*,end=113,err=113) nterm(i,j),nlor(i,j)
 404           v0ij=0.0d0
 405           si=-1.0d0
 406           do k=1,nterm(i,j)
 407             read (itorp,*,end=113,err=113) kk,v1(k,i,j),v2(k,i,j)
 408             v0ij=v0ij+si*v1(k,i,j)
 409             si=-si
 410           enddo
 411           do k=1,nlor(i,j)
 412             read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),
 413      &        vlor2(k,i,j),vlor3(k,i,j)
 414             v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
 415           enddo
 416           v0(i,j)=v0ij
 417         enddo
 418       enddo
 419       close (itorp)
 420       if (lprint) then
 421         write (iout,'(/a/)') 'Torsional constants:'
 422         do i=1,ntortyp
 423           do j=1,ntortyp
 424             write (iout,*) 'ityp',i,' jtyp',j
 425             write (iout,*) 'Fourier constants'
 426             do k=1,nterm(i,j)
 427               write (iout,'(2(1pe15.5))') v1(k,i,j),v2(k,i,j)
 428             enddo
 429             write (iout,*) 'Lorenz constants'
 430             do k=1,nlor(i,j)
 431               write (iout,'(3(1pe15.5))')
 432      &         vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
 433             enddo
 434           enddo
 435         enddo
 436       endif
 437 C
 438 C 6/23/01 Read parameters for double torsionals
 439 C
 440       do i=1,ntortyp
 441         do j=1,ntortyp
 442           do k=1,ntortyp
 443             read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
 444             if (t1.ne.onelett(i) .or. t2.ne.onelett(j)
 445      &        .or. t3.ne.onelett(k)) then
 446               write (iout,*) "Error in double torsional parameter file",
 447      &         i,j,k,t1,t2,t3
 448 #ifdef MPI
 449               call MPI_Finalize(Ierror)
 450 #endif
 451                stop "Error in double torsional parameter file"
 452             endif
 453             read (itordp,*,end=114,err=114) ntermd_1(i,j,k),
 454      &         ntermd_2(i,j,k)
 455             read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k),l=1,
 456      &         ntermd_1(i,j,k))
 457             read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k),l=1,
 458      &         ntermd_1(i,j,k))
 459             read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k),l=1,
 460      &         ntermd_1(i,j,k))
 461             read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k),l=1,
 462      &         ntermd_1(i,j,k))
 463             read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k),
 464      &         v2c(m,l,i,j,k),v2s(l,m,i,j,k),v2s(m,l,i,j,k),
 465      &         m=1,l-1),l=1,ntermd_2(i,j,k))
 466           enddo
 467         enddo
 468       enddo
 469       if (lprint) then
 470       write (iout,*)
 471       write (iout,*) 'Constants for double torsionals'
 472       do i=1,ntortyp
 473         do j=1,ntortyp
 474           do k=1,ntortyp
 475             write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
 476      &        ' nsingle',ntermd_1(i,j,k),' ndouble',ntermd_2(i,j,k)
 477             write (iout,*)
 478             write (iout,*) 'Single angles:'
 479             do l=1,ntermd_1(i,j,k)
 480               write (iout,'(i5,2f10.5,5x,2f10.5)') l,
 481      &           v1c(1,l,i,j,k),v1s(1,l,i,j,k),
 482      &           v1c(2,l,i,j,k),v1s(2,l,i,j,k)
 483             enddo
 484             write (iout,*)
 485             write (iout,*) 'Pairs of angles:'
 486             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k))
 487             do l=1,ntermd_2(i,j,k)
 488               write (iout,'(i5,20f10.5)')
 489      &         l,(v2c(l,m,i,j,k),m=1,ntermd_2(i,j,k))
 490             enddo
 491             write (iout,*)
 492             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k))
 493             do l=1,ntermd_2(i,j,k)
 494               write (iout,'(i5,20f10.5)')
 495      &         l,(v2s(l,m,i,j,k),m=1,ntermd_2(i,j,k))
 496             enddo
 497             write (iout,*)
 498           enddo
 499         enddo
 500       enddo
 501       endif
 502 #endif
 503 C
 504 C 5/21/07 (AL) Read coefficients of the backbone-local sidechain-local
 505 C         correlation energies.
 506 C
 507       read (isccor,*,end=119,err=119) nterm_sccor
 508       do i=1,20
 509         do j=1,20
 510           read (isccor,'(a)')
 511           do k=1,nterm_sccor
 512             read (isccor,*,end=119,err=119) kk,v1sccor(k,i,j),
 513      &        v2sccor(k,i,j)
 514           enddo
 515         enddo
 516       enddo
 517       close (isccor)
 518       if (lprint) then
 519         write (iout,'(/a/)') 'Torsional constants of SCCORR:'
 520         do i=1,20
 521           do j=1,20
 522             write (iout,*) 'ityp',i,' jtyp',j
 523             do k=1,nterm_sccor
 524               write (iout,'(2(1pe15.5))') v1sccor(k,i,j),v2sccor(k,i,j)
 525             enddo
 526           enddo
 527         enddo
 528       endif
 529 C
 530 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
 531 C         interaction energy of the Gly, Ala, and Pro prototypes.
 532 C
 533       if (lprint) then
 534         write (iout,*)
 535         write (iout,*) "Coefficients of the cumulants"
 536       endif
 537       read (ifourier,*) nloctyp
 538       do i=1,nloctyp
 539         read (ifourier,*,end=115,err=115)
 540         read (ifourier,*,end=115,err=115) (b(ii),ii=1,13)
 541         if (lprint) then
 542         write (iout,*) 'Type',i
 543         write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
 544         endif
 545         B1(1,i)  = b(3)
 546         B1(2,i)  = b(5)
 547 c        b1(1,i)=0.0d0
 548 c        b1(2,i)=0.0d0
 549         B1tilde(1,i) = b(3)
 550         B1tilde(2,i) =-b(5)
 551 c        b1tilde(1,i)=0.0d0
 552 c        b1tilde(2,i)=0.0d0
 553         B2(1,i)  = b(2)
 554         B2(2,i)  = b(4)
 555 c        b2(1,i)=0.0d0
 556 c        b2(2,i)=0.0d0
 557         CC(1,1,i)= b(7)
 558         CC(2,2,i)=-b(7)
 559         CC(2,1,i)= b(9)
 560         CC(1,2,i)= b(9)
 561 c        CC(1,1,i)=0.0d0
 562 c        CC(2,2,i)=0.0d0
 563 c        CC(2,1,i)=0.0d0
 564 c        CC(1,2,i)=0.0d0
 565         Ctilde(1,1,i)=b(7)
 566         Ctilde(1,2,i)=b(9)
 567         Ctilde(2,1,i)=-b(9)
 568         Ctilde(2,2,i)=b(7)
 569 c        Ctilde(1,1,i)=0.0d0
 570 c        Ctilde(1,2,i)=0.0d0
 571 c        Ctilde(2,1,i)=0.0d0
 572 c        Ctilde(2,2,i)=0.0d0
 573         DD(1,1,i)= b(6)
 574         DD(2,2,i)=-b(6)
 575         DD(2,1,i)= b(8)
 576         DD(1,2,i)= b(8)
 577 c        DD(1,1,i)=0.0d0
 578 c        DD(2,2,i)=0.0d0
 579 c        DD(2,1,i)=0.0d0
 580 c        DD(1,2,i)=0.0d0
 581         Dtilde(1,1,i)=b(6)
 582         Dtilde(1,2,i)=b(8)
 583         Dtilde(2,1,i)=-b(8)
 584         Dtilde(2,2,i)=b(6)
 585 c        Dtilde(1,1,i)=0.0d0
 586 c        Dtilde(1,2,i)=0.0d0
 587 c        Dtilde(2,1,i)=0.0d0
 588 c        Dtilde(2,2,i)=0.0d0
 589         EE(1,1,i)= b(10)+b(11)
 590         EE(2,2,i)=-b(10)+b(11)
 591         EE(2,1,i)= b(12)-b(13)
 592         EE(1,2,i)= b(12)+b(13)
 593 c        ee(1,1,i)=1.0d0
 594 c        ee(2,2,i)=1.0d0
 595 c        ee(2,1,i)=0.0d0
 596 c        ee(1,2,i)=0.0d0
 597 c        ee(2,1,i)=ee(1,2,i)
 598       enddo
 599       if (lprint) then
 600       do i=1,nloctyp
 601         write (iout,*) 'Type',i
 602         write (iout,*) 'B1'
 603         write(iout,*) B1(1,i),B1(2,i)
 604         write (iout,*) 'B2'
 605         write(iout,*) B2(1,i),B2(2,i)
 606         write (iout,*) 'CC'
 607         do j=1,2
 608           write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
 609         enddo
 610         write(iout,*) 'DD'
 611         do j=1,2
 612           write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
 613         enddo
 614         write(iout,*) 'EE'
 615         do j=1,2
 616           write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
 617         enddo
 618       enddo
 619       endif
 620 C
 621 C Read electrostatic-interaction parameters
 622 C
 623       if (lprint) then
 624         write (iout,*)
 625         write (iout,'(/a)') 'Electrostatic interaction constants:'
 626         write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
 627      &            'IT','JT','APP','BPP','AEL6','AEL3'
 628       endif
 629       read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
 630       read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
 631       read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
 632       read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
 633       close (ielep)
 634       do i=1,2
 635         do j=1,2
 636         rri=rpp(i,j)**6
 637         app (i,j)=epp(i,j)*rri*rri
 638         bpp (i,j)=-2.0D0*epp(i,j)*rri
 639         ael6(i,j)=elpp6(i,j)*4.2D0**6
 640         ael3(i,j)=elpp3(i,j)*4.2D0**3
 641         if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
 642      &                    ael6(i,j),ael3(i,j)
 643         enddo
 644       enddo
 645 C
 646 C Read side-chain interaction parameters.
 647 C
 648       read (isidep,*,end=117,err=117) ipot,expon
 649       if (ipot.lt.1 .or. ipot.gt.5) then
 650         write (iout,'(2a)') 'Error while reading SC interaction',
 651      &               'potential file - unknown potential type.'
 652 #ifdef MPI
 653         call MPI_Finalize(Ierror)
 654 #endif
 655         stop
 656       endif
 657       expon2=expon/2
 658       if(me.eq.king)
 659      & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
 660      & ', exponents are ',expon,2*expon
 661       goto (10,20,30,30,40) ipot
 662 C----------------------- LJ potential ---------------------------------
 663    10 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
 664      &   (sigma0(i),i=1,ntyp)
 665       if (lprint) then
 666         write (iout,'(/a/)') 'Parameters of the LJ potential:'
 667         write (iout,'(a/)') 'The epsilon array:'
 668         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 669         write (iout,'(/a)') 'One-body parameters:'
 670         write (iout,'(a,4x,a)') 'residue','sigma'
 671         write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
 672       endif
 673       goto 50
 674 C----------------------- LJK potential --------------------------------
 675    20 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
 676      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
 677       if (lprint) then
 678         write (iout,'(/a/)') 'Parameters of the LJK potential:'
 679         write (iout,'(a/)') 'The epsilon array:'
 680         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 681         write (iout,'(/a)') 'One-body parameters:'
 682         write (iout,'(a,4x,2a)') 'residue','   sigma  ','    r0    '
 683         write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
 684      &        i=1,ntyp)
 685       endif
 686       goto 50
 687 C---------------------- GB or BP potential -----------------------------
 688    30 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
 689      &  (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip(i),i=1,ntyp),
 690      &  (alp(i),i=1,ntyp)
 691 C For the GB potential convert sigma'**2 into chi'
 692       if (ipot.eq.4) then
 693         do i=1,ntyp
 694           chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
 695         enddo
 696       endif
 697       if (lprint) then
 698         write (iout,'(/a/)') 'Parameters of the BP potential:'
 699         write (iout,'(a/)') 'The epsilon array:'
 700         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 701         write (iout,'(/a)') 'One-body parameters:'
 702         write (iout,'(a,4x,4a)') 'residue','   sigma  ','s||/s_|_^2',
 703      &       '    chip  ','    alph  '
 704         write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
 705      &                     chip(i),alp(i),i=1,ntyp)
 706       endif
 707       goto 50
 708 C--------------------- GBV potential -----------------------------------
 709    40 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
 710      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
 711      &  (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
 712       if (lprint) then
 713         write (iout,'(/a/)') 'Parameters of the GBV potential:'
 714         write (iout,'(a/)') 'The epsilon array:'
 715         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 716         write (iout,'(/a)') 'One-body parameters:'
 717         write (iout,'(a,4x,5a)') 'residue','   sigma  ','    r0    ',
 718      &      's||/s_|_^2','    chip  ','    alph  '
 719         write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
 720      &           sigii(i),chip(i),alp(i),i=1,ntyp)
 721       endif
 722    50 continue
 723       close (isidep)
 724 C-----------------------------------------------------------------------
 725 C Calculate the "working" parameters of SC interactions.
 726       do i=2,ntyp
 727         do j=1,i-1
 728           eps(i,j)=eps(j,i)
 729         enddo
 730       enddo
 731       do i=1,ntyp
 732         do j=i,ntyp
 733           sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
 734           sigma(j,i)=sigma(i,j)
 735           rs0(i,j)=dwa16*sigma(i,j)
 736           rs0(j,i)=rs0(i,j)
 737         enddo
 738       enddo
 739       if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
 740      & 'Working parameters of the SC interactions:',
 741      & '     a    ','     b    ','   augm   ','  sigma ','   r0   ',
 742      & '  chi1   ','   chi2   '
 743       do i=1,ntyp
 744         do j=i,ntyp
 745           epsij=eps(i,j)
 746           if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
 747             rrij=sigma(i,j)
 748           else
 749             rrij=rr0(i)+rr0(j)
 750           endif
 751           r0(i,j)=rrij
 752           r0(j,i)=rrij
 753           rrij=rrij**expon
 754           epsij=eps(i,j)
 755           sigeps=dsign(1.0D0,epsij)
 756           epsij=dabs(epsij)
 757           aa(i,j)=epsij*rrij*rrij
 758           bb(i,j)=-sigeps*epsij*rrij
 759           aa(j,i)=aa(i,j)
 760           bb(j,i)=bb(i,j)
 761           if (ipot.gt.2) then
 762             sigt1sq=sigma0(i)**2
 763             sigt2sq=sigma0(j)**2
 764             sigii1=sigii(i)
 765             sigii2=sigii(j)
 766             ratsig1=sigt2sq/sigt1sq
 767             ratsig2=1.0D0/ratsig1
 768             chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
 769             if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
 770             rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
 771           else
 772             rsum_max=sigma(i,j)
 773           endif
 774 c         if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
 775             sigmaii(i,j)=rsum_max
 776             sigmaii(j,i)=rsum_max
 777 c         else
 778 c           sigmaii(i,j)=r0(i,j)
 779 c           sigmaii(j,i)=r0(i,j)
 780 c         endif
 781 cd        write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
 782           if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
 783             r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
 784             augm(i,j)=epsij*r_augm**(2*expon)
 785 c           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
 786             augm(j,i)=augm(i,j)
 787           else
 788             augm(i,j)=0.0D0
 789             augm(j,i)=0.0D0
 790           endif
 791           if (lprint) then
 792             write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
 793      &      restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
 794      &      sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
 795           endif
 796         enddo
 797       enddo
 798 #ifdef OLDSCP
 799 C
 800 C Define the SC-p interaction constants (hard-coded; old style)
 801 C
 802       do i=1,20
 803 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
 804 C helix formation)
 805 c       aad(i,1)=0.3D0*4.0D0**12
 806 C Following line for constants currently implemented
 807 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
 808         aad(i,1)=1.5D0*4.0D0**12
 809 c       aad(i,1)=0.17D0*5.6D0**12
 810         aad(i,2)=aad(i,1)
 811 C "Soft" SC-p repulsion
 812         bad(i,1)=0.0D0
 813 C Following line for constants currently implemented
 814 c       aad(i,1)=0.3D0*4.0D0**6
 815 C "Hard" SC-p repulsion
 816         bad(i,1)=3.0D0*4.0D0**6
 817 c       bad(i,1)=-2.0D0*0.17D0*5.6D0**6
 818         bad(i,2)=bad(i,1)
 819 c       aad(i,1)=0.0D0
 820 c       aad(i,2)=0.0D0
 821 c       bad(i,1)=1228.8D0
 822 c       bad(i,2)=1228.8D0
 823       enddo
 824 #else
 825 C
 826 C 8/9/01 Read the SC-p interaction constants from file
 827 C
 828       do i=1,ntyp
 829         read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
 830       enddo
 831       do i=1,ntyp
 832         aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
 833         aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
 834         bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
 835         bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
 836       enddo
 837
 838       if (lprint) then
 839         write (iout,*) "Parameters of SC-p interactions:"
 840         do i=1,20
 841           write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
 842      &     eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
 843         enddo
 844       endif
 845 #endif
 846 C
 847 C Define the constants of the disulfide bridge
 848 C
 849       ebr=-5.50D0
 850 c
 851 c Old arbitrary potential - commented out.
 852 c
 853 c      dbr= 4.20D0
 854 c      fbr= 3.30D0
 855 c
 856 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
 857 c energy surface of diethyl disulfide.
 858 c A. Liwo and U. Kozlowska, 11/24/03
 859 c
 860       D0CM = 3.78d0
 861       AKCM = 15.1d0
 862       AKTH = 11.0d0
 863       AKCT = 12.0d0
 864       V1SS =-1.08d0
 865       V2SS = 7.61d0
 866       V3SS = 13.7d0
 867 c      akcm=0.0d0
 868 c      akth=0.0d0
 869 c      akct=0.0d0
 870 c      v1ss=0.0d0
 871 c      v2ss=0.0d0
 872 c      v3ss=0.0d0
 873
 874       if(me.eq.king) then
 875       write (iout,'(/a)') "Disulfide bridge parameters:"
 876       write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
 877       write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
 878       write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
 879       write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
 880      &  ' v3ss:',v3ss
 881       endif
 882       return
 883   111 write (iout,*) "Error reading bending energy parameters."
 884       goto 999
 885   112 write (iout,*) "Error reading rotamer energy parameters."
 886       goto 999
 887   113 write (iout,*) "Error reading torsional energy parameters."
 888       goto 999
 889   114 write (iout,*) "Error reading double torsional energy parameters."
 890       goto 999
 891   115 write (iout,*)
 892      &  "Error reading cumulant (multibody energy) parameters."
 893       goto 999
 894   116 write (iout,*) "Error reading electrostatic energy parameters."
 895       goto 999
 896   117 write (iout,*) "Error reading side chain interaction parameters."
 897       goto 999
 898   118 write (iout,*) "Error reading SCp interaction parameters."
 899       goto 999
 900   119 write (iout,*) "Error reading SCCOR parameters"
 901   999 continue
 902 #ifdef MPI
 903       call MPI_Finalize(Ierror)
 904 #endif
 905       stop
 906       return
 907       end
 908
 909
 910       subroutine getenv_loc(var, val)
 911       character(*) var, val
 912
 913 #ifdef WINIFL
 914       character(2000) line
 915       external ilen
 916
 917       open (196,file='env',status='old',readonly,shared)
 918       iread=0
 919 c      write(*,*)'looking for ',var
 920 10    read(196,*,err=11,end=11)line
 921       iread=index(line,var)
 922 c      write(*,*)iread,' ',var,' ',line
 923       if (iread.eq.0) go to 10
 924 c      write(*,*)'---> ',line
 925 11    continue
 926       if(iread.eq.0) then
 927 c       write(*,*)'CHUJ'
 928        val=''
 929       else
 930        iread=iread+ilen(var)+1
 931        read (line(iread:),*,err=12,end=12) val
 932 c       write(*,*)'OK: ',var,' = ',val
 933       endif
 934       close(196)
 935       return
 936 12    val=''
 937       close(196)
 938 #elif (defined CRAY)
 939       integer lennam,lenval,ierror
 940 c
 941 c        getenv using a POSIX call, useful on the T3D
 942 c        Sept 1996, comment out error check on advice of H. Pritchard
 943 c
 944       lennam = len(var)
 945       if(lennam.le.0) stop '--error calling getenv--'
 946       call pxfgetenv(var,lennam,val,lenval,ierror)
 947 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
 948 #else
 949       call getenv(var,val)
 950 #endif
 951
 952       return
 953       end