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