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       write (2,*) "Start reading THETA_PDB"
 280       do i=1,ntyp
 281         read (ithep_pdb,*,err=111,end=111) a0thet(i),(athet(j,i),j=1,2),
 282      &    (bthet(j,i),j=1,2)
 283         read (ithep_pdb,*,err=111,end=111) (polthet(j,i),j=0,3)
 284         read (ithep_pdb,*,err=111,end=111) (gthet(j,i),j=1,3)
 285         read (ithep_pdb,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
 286         sigc0(i)=sigc0(i)**2
 287       enddo
 288       write (2,*) "End reading THETA_PDB"
 289       close (ithep_pdb)
 290 #endif
 291       close(ithep)
 292 #ifdef CRYST_SC
 293 C
 294 C Read the parameters of the probability distribution/energy expression
 295 C of the side chains.
 296 C
 297       do i=1,ntyp
 298         read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
 299         if (i.eq.10) then
 300           dsc_inv(i)=0.0D0
 301         else
 302           dsc_inv(i)=1.0D0/dsc(i)
 303         endif
 304         if (i.ne.10) then
 305         do j=1,nlob(i)
 306           do k=1,3
 307             do l=1,3
 308               blower(l,k,j)=0.0D0
 309             enddo
 310           enddo
 311         enddo
 312         bsc(1,i)=0.0D0
 313         read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),
 314      &    ((blower(k,l,1),l=1,k),k=1,3)
 315         do j=2,nlob(i)
 316           read (irotam,*,end=112,err=112) bsc(j,i)
 317           read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),
 318      &                                 ((blower(k,l,j),l=1,k),k=1,3)
 319         enddo
 320         do j=1,nlob(i)
 321           do k=1,3
 322             do l=1,k
 323               akl=0.0D0
 324               do m=1,3
 325                 akl=akl+blower(k,m,j)*blower(l,m,j)
 326               enddo
 327               gaussc(k,l,j,i)=akl
 328               gaussc(l,k,j,i)=akl
 329             enddo
 330           enddo
 331         enddo
 332         endif
 333       enddo
 334       close (irotam)
 335       if (lprint) then
 336         write (iout,'(/a)') 'Parameters of side-chain local geometry'
 337         do i=1,ntyp
 338           nlobi=nlob(i)
 339           if (nlobi.gt.0) then
 340             if (LaTeX) then
 341               write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
 342      &         ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
 343                write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
 344      &                             'log h',(bsc(j,i),j=1,nlobi)
 345                write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
 346      &        'x',((censc(k,j,i),k=1,3),j=1,nlobi)
 347               do k=1,3
 348                 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
 349      &                 ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
 350               enddo
 351             else
 352               write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
 353               write (iout,'(a,f10.4,4(16x,f10.4))')
 354      &                             'Center  ',(bsc(j,i),j=1,nlobi)
 355               write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),
 356      &           j=1,nlobi)
 357               write (iout,'(a)')
 358             endif
 359           endif
 360         enddo
 361       endif
 362 #else
 363 C
 364 C Read scrot parameters for potentials determined from all-atom AM1 calculations
 365 C added by Urszula Kozlowska 07/11/2007
 366 C
 367       do i=1,ntyp
 368         read (irotam,*,end=112,err=112)
 369        if (i.eq.10) then
 370          read (irotam,*,end=112,err=112)
 371        else
 372          do j=1,65
 373            read(irotam,*,end=112,err=112) sc_parmin(j,i)
 374          enddo
 375        endif
 376       enddo
 377 C
 378 C Read the parameters of the probability distribution/energy expression
 379 C of the side chains.
 380 C
 381       do i=1,ntyp
 382         read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
 383         if (i.eq.10) then
 384           dsc_inv(i)=0.0D0
 385         else
 386           dsc_inv(i)=1.0D0/dsc(i)
 387         endif
 388         if (i.ne.10) then
 389         do j=1,nlob(i)
 390           do k=1,3
 391             do l=1,3
 392               blower(l,k,j)=0.0D0
 393             enddo
 394           enddo
 395         enddo
 396         bsc(1,i)=0.0D0
 397         read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),
 398      &    ((blower(k,l,1),l=1,k),k=1,3)
 399         do j=2,nlob(i)
 400           read (irotam_pdb,*,end=112,err=112) bsc(j,i)
 401           read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),
 402      &                                 ((blower(k,l,j),l=1,k),k=1,3)
 403         enddo
 404         do j=1,nlob(i)
 405           do k=1,3
 406             do l=1,k
 407               akl=0.0D0
 408               do m=1,3
 409                 akl=akl+blower(k,m,j)*blower(l,m,j)
 410               enddo
 411               gaussc(k,l,j,i)=akl
 412               gaussc(l,k,j,i)=akl
 413             enddo
 414           enddo
 415         enddo
 416         endif
 417       enddo
 418       close (irotam_pdb)
 419 #endif
 420       close(irotam)
 421
 422 #ifdef CRYST_TOR
 423 C
 424 C Read torsional parameters in old format
 425 C
 426       read (itorp,*,end=113,err=113) ntortyp,nterm_old
 427       if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
 428       read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
 429       do i=1,ntortyp
 430         do j=1,ntortyp
 431           read (itorp,'(a)')
 432           do k=1,nterm_old
 433             read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
 434           enddo
 435         enddo
 436       enddo
 437       close (itorp)
 438       if (lprint) then
 439         write (iout,'(/a/)') 'Torsional constants:'
 440         do i=1,ntortyp
 441           do j=1,ntortyp
 442             write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
 443             write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
 444           enddo
 445         enddo
 446       endif
 447 #else
 448 C
 449 C Read torsional parameters
 450 C
 451       read (itorp,*,end=113,err=113) ntortyp
 452       read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
 453 c      write (iout,*) 'ntortyp',ntortyp
 454       do i=1,ntortyp
 455         do j=1,ntortyp
 456           read (itorp,*,end=113,err=113) nterm(i,j),nlor(i,j)
 457           v0ij=0.0d0
 458           si=-1.0d0
 459           do k=1,nterm(i,j)
 460             read (itorp,*,end=113,err=113) kk,v1(k,i,j),v2(k,i,j)
 461             v0ij=v0ij+si*v1(k,i,j)
 462             si=-si
 463           enddo
 464           do k=1,nlor(i,j)
 465             read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),
 466      &        vlor2(k,i,j),vlor3(k,i,j)
 467             v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
 468           enddo
 469           v0(i,j)=v0ij
 470         enddo
 471       enddo
 472       close (itorp)
 473       if (lprint) then
 474         write (iout,'(/a/)') 'Torsional constants:'
 475         do i=1,ntortyp
 476           do j=1,ntortyp
 477             write (iout,*) 'ityp',i,' jtyp',j
 478             write (iout,*) 'Fourier constants'
 479             do k=1,nterm(i,j)
 480               write (iout,'(2(1pe15.5))') v1(k,i,j),v2(k,i,j)
 481             enddo
 482             write (iout,*) 'Lorenz constants'
 483             do k=1,nlor(i,j)
 484               write (iout,'(3(1pe15.5))')
 485      &         vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
 486             enddo
 487           enddo
 488         enddo
 489       endif
 490 C
 491 C 6/23/01 Read parameters for double torsionals
 492 C
 493       do i=1,ntortyp
 494         do j=1,ntortyp
 495           do k=1,ntortyp
 496             read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
 497             if (t1.ne.onelett(i) .or. t2.ne.onelett(j)
 498      &        .or. t3.ne.onelett(k)) then
 499               write (iout,*) "Error in double torsional parameter file",
 500      &         i,j,k,t1,t2,t3
 501 #ifdef MPI
 502               call MPI_Finalize(Ierror)
 503 #endif
 504                stop "Error in double torsional parameter file"
 505             endif
 506             read (itordp,*,end=114,err=114) ntermd_1(i,j,k),
 507      &         ntermd_2(i,j,k)
 508             read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k),l=1,
 509      &         ntermd_1(i,j,k))
 510             read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k),l=1,
 511      &         ntermd_1(i,j,k))
 512             read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k),l=1,
 513      &         ntermd_1(i,j,k))
 514             read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k),l=1,
 515      &         ntermd_1(i,j,k))
 516             read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k),
 517      &         v2c(m,l,i,j,k),v2s(l,m,i,j,k),v2s(m,l,i,j,k),
 518      &         m=1,l-1),l=1,ntermd_2(i,j,k))
 519           enddo
 520         enddo
 521       enddo
 522       if (lprint) then
 523       write (iout,*)
 524       write (iout,*) 'Constants for double torsionals'
 525       do i=1,ntortyp
 526         do j=1,ntortyp
 527           do k=1,ntortyp
 528             write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
 529      &        ' nsingle',ntermd_1(i,j,k),' ndouble',ntermd_2(i,j,k)
 530             write (iout,*)
 531             write (iout,*) 'Single angles:'
 532             do l=1,ntermd_1(i,j,k)
 533               write (iout,'(i5,2f10.5,5x,2f10.5)') l,
 534      &           v1c(1,l,i,j,k),v1s(1,l,i,j,k),
 535      &           v1c(2,l,i,j,k),v1s(2,l,i,j,k)
 536             enddo
 537             write (iout,*)
 538             write (iout,*) 'Pairs of angles:'
 539             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k))
 540             do l=1,ntermd_2(i,j,k)
 541               write (iout,'(i5,20f10.5)')
 542      &         l,(v2c(l,m,i,j,k),m=1,ntermd_2(i,j,k))
 543             enddo
 544             write (iout,*)
 545             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k))
 546             do l=1,ntermd_2(i,j,k)
 547               write (iout,'(i5,20f10.5)')
 548      &         l,(v2s(l,m,i,j,k),m=1,ntermd_2(i,j,k))
 549             enddo
 550             write (iout,*)
 551           enddo
 552         enddo
 553       enddo
 554       endif
 555 #endif
 556 C Read of Side-chain backbone correlation parameters
 557 C Modified 11 May 2012 by Adasko
 558 CCC
 559 C
 560       read (isccor,*,end=113,err=113) nsccortyp
 561       read (isccor,*,end=113,err=113) (isccortyp(i),i=1,ntyp)
 562 c      write (iout,*) 'ntortyp',ntortyp
 563       maxinter=3
 564 cc maxinter is maximum interaction sites
 565       do l=1,maxinter
 566       do i=1,nsccortyp
 567         do j=1,nsccortyp
 568           read (isccor,*,end=113,err=113) nterm_sccor(i,j),nlor_sccor(i,j)
 569           v0ijsccor=0.0d0
 570           si=-1.0d0
 571           do k=1,nterm_sccor(i,j)
 572             read (isccor,*,end=113,err=113) kk,v1sccor(k,i,j),v2sccor(k,i,j)
 573             v0ijsccor=v0ijsccor+si*v1sccor(k,i,j)
 574             si=-si
 575           enddo
 576           do k=1,nlor_sccor(i,j)
 577             read (isccor,*,end=113,err=113) kk,vlor1sccor(k,i,j),
 578      &        vlor2sccor(k,i,j),vlor3sccor(k,i,j)
 579             v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
 580      &(1+vlor3sccor(k,i,j)**2)
 581           enddo
 582           v0sccor(i,j)=v0ijsccor
 583         enddo
 584       enddo
 585       enddo
 586       close (isccor)
 587
 588       if (lprint) then
 589         write (iout,'(/a/)') 'Torsional constants:'
 590         do i=1,nsccortyp
 591           do j=1,nsccortyp
 592             write (iout,*) 'ityp',i,' jtyp',j
 593             write (iout,*) 'Fourier constants'
 594             do k=1,nterm_sccor(i,j)
 595               write (iout,'(2(1pe15.5))') v1sccor(k,i,j),v2sccor(k,i,j)
 596             enddo
 597             write (iout,*) 'Lorenz constants'
 598             do k=1,nlor_sccor(i,j)
 599               write (iout,'(3(1pe15.5))')
 600      &         vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
 601             enddo
 602           enddo
 603         enddo
 604       endif
 605 C
 606 C
 607 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
 608 C         interaction energy of the Gly, Ala, and Pro prototypes.
 609 C
 610       if (lprint) then
 611         write (iout,*)
 612         write (iout,*) "Coefficients of the cumulants"
 613       endif
 614       read (ifourier,*) nloctyp
 615       do i=1,nloctyp
 616         read (ifourier,*,end=115,err=115)
 617         read (ifourier,*,end=115,err=115) (b(ii),ii=1,13)
 618         if (lprint) then
 619         write (iout,*) 'Type',i
 620         write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
 621         endif
 622         B1(1,i)  = b(3)
 623         B1(2,i)  = b(5)
 624 c        b1(1,i)=0.0d0
 625 c        b1(2,i)=0.0d0
 626         B1tilde(1,i) = b(3)
 627         B1tilde(2,i) =-b(5)
 628 c        b1tilde(1,i)=0.0d0
 629 c        b1tilde(2,i)=0.0d0
 630         B2(1,i)  = b(2)
 631         B2(2,i)  = b(4)
 632 c        b2(1,i)=0.0d0
 633 c        b2(2,i)=0.0d0
 634         CC(1,1,i)= b(7)
 635         CC(2,2,i)=-b(7)
 636         CC(2,1,i)= b(9)
 637         CC(1,2,i)= b(9)
 638 c        CC(1,1,i)=0.0d0
 639 c        CC(2,2,i)=0.0d0
 640 c        CC(2,1,i)=0.0d0
 641 c        CC(1,2,i)=0.0d0
 642         Ctilde(1,1,i)=b(7)
 643         Ctilde(1,2,i)=b(9)
 644         Ctilde(2,1,i)=-b(9)
 645         Ctilde(2,2,i)=b(7)
 646 c        Ctilde(1,1,i)=0.0d0
 647 c        Ctilde(1,2,i)=0.0d0
 648 c        Ctilde(2,1,i)=0.0d0
 649 c        Ctilde(2,2,i)=0.0d0
 650         DD(1,1,i)= b(6)
 651         DD(2,2,i)=-b(6)
 652         DD(2,1,i)= b(8)
 653         DD(1,2,i)= b(8)
 654 c        DD(1,1,i)=0.0d0
 655 c        DD(2,2,i)=0.0d0
 656 c        DD(2,1,i)=0.0d0
 657 c        DD(1,2,i)=0.0d0
 658         Dtilde(1,1,i)=b(6)
 659         Dtilde(1,2,i)=b(8)
 660         Dtilde(2,1,i)=-b(8)
 661         Dtilde(2,2,i)=b(6)
 662 c        Dtilde(1,1,i)=0.0d0
 663 c        Dtilde(1,2,i)=0.0d0
 664 c        Dtilde(2,1,i)=0.0d0
 665 c        Dtilde(2,2,i)=0.0d0
 666         EE(1,1,i)= b(10)+b(11)
 667         EE(2,2,i)=-b(10)+b(11)
 668         EE(2,1,i)= b(12)-b(13)
 669         EE(1,2,i)= b(12)+b(13)
 670 c        ee(1,1,i)=1.0d0
 671 c        ee(2,2,i)=1.0d0
 672 c        ee(2,1,i)=0.0d0
 673 c        ee(1,2,i)=0.0d0
 674 c        ee(2,1,i)=ee(1,2,i)
 675       enddo
 676       if (lprint) then
 677       do i=1,nloctyp
 678         write (iout,*) 'Type',i
 679         write (iout,*) 'B1'
 680         write(iout,*) B1(1,i),B1(2,i)
 681         write (iout,*) 'B2'
 682         write(iout,*) B2(1,i),B2(2,i)
 683         write (iout,*) 'CC'
 684         do j=1,2
 685           write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
 686         enddo
 687         write(iout,*) 'DD'
 688         do j=1,2
 689           write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
 690         enddo
 691         write(iout,*) 'EE'
 692         do j=1,2
 693           write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
 694         enddo
 695       enddo
 696       endif
 697 C
 698 C Read electrostatic-interaction parameters
 699 C
 700       if (lprint) then
 701         write (iout,*)
 702         write (iout,'(/a)') 'Electrostatic interaction constants:'
 703         write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
 704      &            'IT','JT','APP','BPP','AEL6','AEL3'
 705       endif
 706       read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
 707       read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
 708       read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
 709       read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
 710       close (ielep)
 711       do i=1,2
 712         do j=1,2
 713         rri=rpp(i,j)**6
 714         app (i,j)=epp(i,j)*rri*rri
 715         bpp (i,j)=-2.0D0*epp(i,j)*rri
 716         ael6(i,j)=elpp6(i,j)*4.2D0**6
 717         ael3(i,j)=elpp3(i,j)*4.2D0**3
 718         if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
 719      &                    ael6(i,j),ael3(i,j)
 720         enddo
 721       enddo
 722 C
 723 C Read side-chain interaction parameters.
 724 C
 725       read (isidep,*,end=117,err=117) ipot,expon
 726       if (ipot.lt.1 .or. ipot.gt.5) then
 727         write (iout,'(2a)') 'Error while reading SC interaction',
 728      &               'potential file - unknown potential type.'
 729 #ifdef MPI
 730         call MPI_Finalize(Ierror)
 731 #endif
 732         stop
 733       endif
 734       expon2=expon/2
 735       if(me.eq.king)
 736      & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
 737      & ', exponents are ',expon,2*expon
 738       goto (10,20,30,30,40) ipot
 739 C----------------------- LJ potential ---------------------------------
 740    10 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
 741      &   (sigma0(i),i=1,ntyp)
 742       if (lprint) then
 743         write (iout,'(/a/)') 'Parameters of the LJ potential:'
 744         write (iout,'(a/)') 'The epsilon array:'
 745         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 746         write (iout,'(/a)') 'One-body parameters:'
 747         write (iout,'(a,4x,a)') 'residue','sigma'
 748         write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
 749       endif
 750       goto 50
 751 C----------------------- LJK potential --------------------------------
 752    20 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
 753      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
 754       if (lprint) then
 755         write (iout,'(/a/)') 'Parameters of the LJK potential:'
 756         write (iout,'(a/)') 'The epsilon array:'
 757         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 758         write (iout,'(/a)') 'One-body parameters:'
 759         write (iout,'(a,4x,2a)') 'residue','   sigma  ','    r0    '
 760         write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
 761      &        i=1,ntyp)
 762       endif
 763       goto 50
 764 C---------------------- GB or BP potential -----------------------------
 765    30 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
 766      &  (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip(i),i=1,ntyp),
 767      &  (alp(i),i=1,ntyp)
 768 C For the GB potential convert sigma'**2 into chi'
 769       if (ipot.eq.4) then
 770         do i=1,ntyp
 771           chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
 772         enddo
 773       endif
 774       if (lprint) then
 775         write (iout,'(/a/)') 'Parameters of the BP potential:'
 776         write (iout,'(a/)') 'The epsilon array:'
 777         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 778         write (iout,'(/a)') 'One-body parameters:'
 779         write (iout,'(a,4x,4a)') 'residue','   sigma  ','s||/s_|_^2',
 780      &       '    chip  ','    alph  '
 781         write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
 782      &                     chip(i),alp(i),i=1,ntyp)
 783       endif
 784       goto 50
 785 C--------------------- GBV potential -----------------------------------
 786    40 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
 787      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
 788      &  (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
 789       if (lprint) then
 790         write (iout,'(/a/)') 'Parameters of the GBV potential:'
 791         write (iout,'(a/)') 'The epsilon array:'
 792         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 793         write (iout,'(/a)') 'One-body parameters:'
 794         write (iout,'(a,4x,5a)') 'residue','   sigma  ','    r0    ',
 795      &      's||/s_|_^2','    chip  ','    alph  '
 796         write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
 797      &           sigii(i),chip(i),alp(i),i=1,ntyp)
 798       endif
 799    50 continue
 800       close (isidep)
 801 C-----------------------------------------------------------------------
 802 C Calculate the "working" parameters of SC interactions.
 803       do i=2,ntyp
 804         do j=1,i-1
 805           eps(i,j)=eps(j,i)
 806         enddo
 807       enddo
 808       do i=1,ntyp
 809         do j=i,ntyp
 810           sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
 811           sigma(j,i)=sigma(i,j)
 812           rs0(i,j)=dwa16*sigma(i,j)
 813           rs0(j,i)=rs0(i,j)
 814         enddo
 815       enddo
 816       if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
 817      & 'Working parameters of the SC interactions:',
 818      & '     a    ','     b    ','   augm   ','  sigma ','   r0   ',
 819      & '  chi1   ','   chi2   '
 820       do i=1,ntyp
 821         do j=i,ntyp
 822           epsij=eps(i,j)
 823           if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
 824             rrij=sigma(i,j)
 825           else
 826             rrij=rr0(i)+rr0(j)
 827           endif
 828           r0(i,j)=rrij
 829           r0(j,i)=rrij
 830           rrij=rrij**expon
 831           epsij=eps(i,j)
 832           sigeps=dsign(1.0D0,epsij)
 833           epsij=dabs(epsij)
 834           aa(i,j)=epsij*rrij*rrij
 835           bb(i,j)=-sigeps*epsij*rrij
 836           aa(j,i)=aa(i,j)
 837           bb(j,i)=bb(i,j)
 838           if (ipot.gt.2) then
 839             sigt1sq=sigma0(i)**2
 840             sigt2sq=sigma0(j)**2
 841             sigii1=sigii(i)
 842             sigii2=sigii(j)
 843             ratsig1=sigt2sq/sigt1sq
 844             ratsig2=1.0D0/ratsig1
 845             chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
 846             if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
 847             rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
 848           else
 849             rsum_max=sigma(i,j)
 850           endif
 851 c         if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
 852             sigmaii(i,j)=rsum_max
 853             sigmaii(j,i)=rsum_max
 854 c         else
 855 c           sigmaii(i,j)=r0(i,j)
 856 c           sigmaii(j,i)=r0(i,j)
 857 c         endif
 858 cd        write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
 859           if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
 860             r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
 861             augm(i,j)=epsij*r_augm**(2*expon)
 862 c           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
 863             augm(j,i)=augm(i,j)
 864           else
 865             augm(i,j)=0.0D0
 866             augm(j,i)=0.0D0
 867           endif
 868           if (lprint) then
 869             write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
 870      &      restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
 871      &      sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
 872           endif
 873         enddo
 874       enddo
 875 #ifdef OLDSCP
 876 C
 877 C Define the SC-p interaction constants (hard-coded; old style)
 878 C
 879       do i=1,20
 880 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
 881 C helix formation)
 882 c       aad(i,1)=0.3D0*4.0D0**12
 883 C Following line for constants currently implemented
 884 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
 885         aad(i,1)=1.5D0*4.0D0**12
 886 c       aad(i,1)=0.17D0*5.6D0**12
 887         aad(i,2)=aad(i,1)
 888 C "Soft" SC-p repulsion
 889         bad(i,1)=0.0D0
 890 C Following line for constants currently implemented
 891 c       aad(i,1)=0.3D0*4.0D0**6
 892 C "Hard" SC-p repulsion
 893         bad(i,1)=3.0D0*4.0D0**6
 894 c       bad(i,1)=-2.0D0*0.17D0*5.6D0**6
 895         bad(i,2)=bad(i,1)
 896 c       aad(i,1)=0.0D0
 897 c       aad(i,2)=0.0D0
 898 c       bad(i,1)=1228.8D0
 899 c       bad(i,2)=1228.8D0
 900       enddo
 901 #else
 902 C
 903 C 8/9/01 Read the SC-p interaction constants from file
 904 C
 905       do i=1,ntyp
 906         read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
 907       enddo
 908       do i=1,ntyp
 909         aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
 910         aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
 911         bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
 912         bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
 913       enddo
 914
 915       if (lprint) then
 916         write (iout,*) "Parameters of SC-p interactions:"
 917         do i=1,20
 918           write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
 919      &     eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
 920         enddo
 921       endif
 922 #endif
 923 C
 924 C Define the constants of the disulfide bridge
 925 C
 926       ebr=-5.50D0
 927 c
 928 c Old arbitrary potential - commented out.
 929 c
 930 c      dbr= 4.20D0
 931 c      fbr= 3.30D0
 932 c
 933 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
 934 c energy surface of diethyl disulfide.
 935 c A. Liwo and U. Kozlowska, 11/24/03
 936 c
 937       D0CM = 3.78d0
 938       AKCM = 15.1d0
 939       AKTH = 11.0d0
 940       AKCT = 12.0d0
 941       V1SS =-1.08d0
 942       V2SS = 7.61d0
 943       V3SS = 13.7d0
 944 c      akcm=0.0d0
 945 c      akth=0.0d0
 946 c      akct=0.0d0
 947 c      v1ss=0.0d0
 948 c      v2ss=0.0d0
 949 c      v3ss=0.0d0
 950
 951       if(me.eq.king) then
 952       write (iout,'(/a)') "Disulfide bridge parameters:"
 953       write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
 954       write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
 955       write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
 956       write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
 957      &  ' v3ss:',v3ss
 958       endif
 959       return
 960   111 write (iout,*) "Error reading bending energy parameters."
 961       goto 999
 962   112 write (iout,*) "Error reading rotamer energy parameters."
 963       goto 999
 964   113 write (iout,*) "Error reading torsional energy parameters."
 965       goto 999
 966   114 write (iout,*) "Error reading double torsional energy parameters."
 967       goto 999
 968   115 write (iout,*)
 969      &  "Error reading cumulant (multibody energy) parameters."
 970       goto 999
 971   116 write (iout,*) "Error reading electrostatic energy parameters."
 972       goto 999
 973   117 write (iout,*) "Error reading side chain interaction parameters."
 974       goto 999
 975   118 write (iout,*) "Error reading SCp interaction parameters."
 976       goto 999
 977   119 write (iout,*) "Error reading SCCOR parameters"
 978   999 continue
 979 #ifdef MPI
 980       call MPI_Finalize(Ierror)
 981 #endif
 982       stop
 983       return
 984       end
 985
 986
 987       subroutine getenv_loc(var, val)
 988       character(*) var, val
 989
 990 #ifdef WINIFL
 991       character(2000) line
 992       external ilen
 993
 994       open (196,file='env',status='old',readonly,shared)
 995       iread=0
 996 c      write(*,*)'looking for ',var
 997 10    read(196,*,err=11,end=11)line
 998       iread=index(line,var)
 999 c      write(*,*)iread,' ',var,' ',line
1000       if (iread.eq.0) go to 10
1001 c      write(*,*)'---> ',line
1002 11    continue
1003       if(iread.eq.0) then
1004 c       write(*,*)'CHUJ'
1005        val=''
1006       else
1007        iread=iread+ilen(var)+1
1008        read (line(iread:),*,err=12,end=12) val
1009 c       write(*,*)'OK: ',var,' = ',val
1010       endif
1011       close(196)
1012       return
1013 12    val=''
1014       close(196)
1015 #elif (defined CRAY)
1016       integer lennam,lenval,ierror
1017 c
1018 c        getenv using a POSIX call, useful on the T3D
1019 c        Sept 1996, comment out error check on advice of H. Pritchard
1020 c
1021       lennam = len(var)
1022       if(lennam.le.0) stop '--error calling getenv--'
1023       call pxfgetenv(var,lennam,val,lenval,ierror)
1024 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1025 #else
1026       call getenv(var,val)
1027 #endif
1028
1029       return
1030       end