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       do iblock=1,2
 455       do i=-ntyp,-1
 456        itortyp(i)=-itortyp(-1)
 457       enddo
 458 c      write (iout,*) 'ntortyp',ntortyp
 459       do i=0,ntortyp,ntortyp-1
 460         do j=-ntortyp,ntortyp
 461           read (itorp,*,end=113,err=113) nterm(i,j,iblock),
 462      &          nlor(i,j,iblock)
 463            nterm(-i,-j,iblock)=nterm(i,j,iblock)
 464            nlor(-i,-j,iblock)=nlor(i,j,iblock)
 465           v0ij=0.0d0
 466           si=-1.0d0
 467           do k=1,nterm(i,j,iblock)
 468             read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),
 469      &      v2(k,i,j,iblock)
 470             v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
 471             v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
 472             v0ij=v0ij+si*v1(k,i,j,iblock)
 473             si=-si
 474           enddo
 475           do k=1,nlor(i,j,iblock)
 476             read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),
 477      &        vlor2(k,i,j),vlor3(k,i,j)
 478             v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
 479           enddo
 480           v0(i,j)=v0ij
 481         enddo
 482       enddo
 483       close (itorp)
 484       if (lprint) then
 485         write (iout,'(/a/)') 'Torsional constants:'
 486         do i=1,ntortyp
 487           do j=1,ntortyp
 488             write (iout,*) 'ityp',i,' jtyp',j
 489             write (iout,*) 'Fourier constants'
 490             do k=1,nterm(i,j,iblock)
 491               write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
 492      &         v2(k,i,j,iblock)
 493             enddo
 494             write (iout,*) 'Lorenz constants'
 495             do k=1,nlor(i,j,iblock)
 496               write (iout,'(3(1pe15.5))')
 497      &         vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
 498             enddo
 499           enddo
 500         enddo
 501       endif
 502 C
 503 C 6/23/01 Read parameters for double torsionals
 504 C
 505       do iblock=1,2
 506       do i=0,ntortyp-1
 507         do j=-ntortyp+1,ntortyp-1
 508           do k=-ntortyp+1,ntortyp-1
 509             read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
 510               write (iout,*) "OK onelett",
 511      &         i,j,k,t1,t2,t3
 512
 513             if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
 514      &        .or. t3.ne.toronelet(k)) then
 515               write (iout,*) "Error in double torsional parameter file",
 516      &         i,j,k,t1,t2,t3
 517 #ifdef MPI
 518               call MPI_Finalize(Ierror)
 519 #endif
 520                stop "Error in double torsional parameter file"
 521             endif
 522             read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),
 523      &         ntermd_2(i,j,k,iblock)
 524             read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,
 525      &         ntermd_1(i,j,k,iblock))
 526             read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,
 527      &         ntermd_1(i,j,k,iblock))
 528             read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,
 529      &         ntermd_1(i,j,k,iblock))
 530             read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,
 531      &         ntermd_1(i,j,k,iblock))
 532 C Matrix of D parameters for one dimesional foureir series
 533             do l=1,  ntermd_1(i,j,k,iblock)
 534              v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
 535              v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
 536              v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
 537              v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
 538             enddo
 539             read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),
 540      &         v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
 541      &         v2s(m,l,i,j,k,iblock),
 542      &         m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
 543 C Matrix of D parameters for two dimesional fourier series
 544             do l=1,ntermd_2(i,j,k,iblock)
 545              do m=1,l-1
 546              v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
 547              v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
 548              v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
 549              v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
 550              enddo!m
 551             enddo!l
 552           enddo!k
 553         enddo!j
 554       enddo!i
 555       enddo!iblock
 556 cc      if (lprint) then
 557       write (iout,*)
 558       write (iout,*) 'Constants for double torsionals'
 559       do i=1,ntortyp
 560         do j=-ntortyp,ntortyp
 561           do k=-ntortyp,ntortyp
 562             write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
 563      &        ' nsingle',ntermd_1(i,j,k,iblock),
 564      &        ' ndouble',ntermd_2(i,j,k,iblock)
 565             write (iout,*)
 566             write (iout,*) 'Single angles:'
 567             do l=1,ntermd_1(i,j,k,iblock)
 568               write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
 569      &           v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
 570      &           v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock)
 571             enddo
 572             write (iout,*)
 573             write (iout,*) 'Pairs of angles:'
 574             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
 575             do l=1,ntermd_2(i,j,k,iblock)
 576               write (iout,'(i5,20f10.5)')
 577      &         l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
 578             enddo
 579             write (iout,*)
 580             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
 581             do l=1,ntermd_2(i,j,k,iblock)
 582               write (iout,'(i5,20f10.5)')
 583      &         l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
 584      &         (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
 585             enddo
 586             write (iout,*)
 587           enddo
 588         enddo
 589       enddo
 590 cc      endif
 591 #endif
 592 C
 593 C 5/21/07 (AL) Read coefficients of the backbone-local sidechain-local
 594 C         correlation energies.
 595 C
 596       read (isccor,*,end=119,err=119) nterm_sccor
 597       do i=1,20
 598         do j=1,20
 599           read (isccor,'(a)')
 600           do k=1,nterm_sccor
 601             read (isccor,*,end=119,err=119) kk,v1sccor(k,i,j),
 602      &        v2sccor(k,i,j)
 603           enddo
 604         enddo
 605       enddo
 606       close (isccor)
 607       if (lprint) then
 608         write (iout,'(/a/)') 'Torsional constants of SCCORR:'
 609         do i=1,20
 610           do j=1,20
 611             write (iout,*) 'ityp',i,' jtyp',j
 612             do k=1,nterm_sccor
 613               write (iout,'(2(1pe15.5))') v1sccor(k,i,j),v2sccor(k,i,j)
 614             enddo
 615           enddo
 616         enddo
 617       endif
 618 C
 619 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
 620 C         interaction energy of the Gly, Ala, and Pro prototypes.
 621 C
 622       if (lprint) then
 623         write (iout,*)
 624         write (iout,*) "Coefficients of the cumulants"
 625       endif
 626       read (ifourier,*) nloctyp
 627       do i=1,nloctyp
 628       do i=0,nloctyp-1
 629         read (ifourier,*,end=115,err=115)
 630         read (ifourier,*,end=115,err=115) (b(ii),ii=1,13)
 631         if (lprint) then
 632         write (iout,*) 'Type',i
 633         write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
 634         endif
 635         B1(1,i)  = b(3)
 636         B1(2,i)  = b(5)
 637 c        b1(1,i)=0.0d0
 638 c        b1(2,i)=0.0d0
 639         B1tilde(1,i) = b(3)
 640         B1tilde(2,i) =-b(5)
 641 c        b1tilde(1,i)=0.0d0
 642 c        b1tilde(2,i)=0.0d0
 643         B2(1,i)  = b(2)
 644         B2(2,i)  = b(4)
 645 c        b2(1,i)=0.0d0
 646 c        b2(2,i)=0.0d0
 647         CC(1,1,i)= b(7)
 648         CC(2,2,i)=-b(7)
 649         CC(2,1,i)= b(9)
 650         CC(1,2,i)= b(9)
 651 c        CC(1,1,i)=0.0d0
 652 c        CC(2,2,i)=0.0d0
 653 c        CC(2,1,i)=0.0d0
 654 c        CC(1,2,i)=0.0d0
 655         Ctilde(1,1,i)=b(7)
 656         Ctilde(1,2,i)=b(9)
 657         Ctilde(2,1,i)=-b(9)
 658         Ctilde(2,2,i)=b(7)
 659 c        Ctilde(1,1,i)=0.0d0
 660 c        Ctilde(1,2,i)=0.0d0
 661 c        Ctilde(2,1,i)=0.0d0
 662 c        Ctilde(2,2,i)=0.0d0
 663         DD(1,1,i)= b(6)
 664         DD(2,2,i)=-b(6)
 665         DD(2,1,i)= b(8)
 666         DD(1,2,i)= b(8)
 667 c        DD(1,1,i)=0.0d0
 668 c        DD(2,2,i)=0.0d0
 669 c        DD(2,1,i)=0.0d0
 670 c        DD(1,2,i)=0.0d0
 671         Dtilde(1,1,i)=b(6)
 672         Dtilde(1,2,i)=b(8)
 673         Dtilde(2,1,i)=-b(8)
 674         Dtilde(2,2,i)=b(6)
 675 c        Dtilde(1,1,i)=0.0d0
 676 c        Dtilde(1,2,i)=0.0d0
 677 c        Dtilde(2,1,i)=0.0d0
 678 c        Dtilde(2,2,i)=0.0d0
 679         EE(1,1,i)= b(10)+b(11)
 680         EE(2,2,i)=-b(10)+b(11)
 681         EE(2,1,i)= b(12)-b(13)
 682         EE(1,2,i)= b(12)+b(13)
 683 c        ee(1,1,i)=1.0d0
 684 c        ee(2,2,i)=1.0d0
 685 c        ee(2,1,i)=0.0d0
 686 c        ee(1,2,i)=0.0d0
 687 c        ee(2,1,i)=ee(1,2,i)
 688       enddo
 689       if (lprint) then
 690       do i=1,nloctyp
 691         write (iout,*) 'Type',i
 692         write (iout,*) 'B1'
 693         write(iout,*) B1(1,i),B1(2,i)
 694         write (iout,*) 'B2'
 695         write(iout,*) B2(1,i),B2(2,i)
 696         write (iout,*) 'CC'
 697         do j=1,2
 698           write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
 699         enddo
 700         write(iout,*) 'DD'
 701         do j=1,2
 702           write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
 703         enddo
 704         write(iout,*) 'EE'
 705         do j=1,2
 706           write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
 707         enddo
 708       enddo
 709       endif
 710 C
 711 C Read electrostatic-interaction parameters
 712 C
 713       if (lprint) then
 714         write (iout,*)
 715         write (iout,'(/a)') 'Electrostatic interaction constants:'
 716         write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
 717      &            'IT','JT','APP','BPP','AEL6','AEL3'
 718       endif
 719       read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
 720       read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
 721       read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
 722       read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
 723       close (ielep)
 724       do i=1,2
 725         do j=1,2
 726         rri=rpp(i,j)**6
 727         app (i,j)=epp(i,j)*rri*rri
 728         bpp (i,j)=-2.0D0*epp(i,j)*rri
 729         ael6(i,j)=elpp6(i,j)*4.2D0**6
 730         ael3(i,j)=elpp3(i,j)*4.2D0**3
 731         if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
 732      &                    ael6(i,j),ael3(i,j)
 733         enddo
 734       enddo
 735 C
 736 C Read side-chain interaction parameters.
 737 C
 738       read (isidep,*,end=117,err=117) ipot,expon
 739       if (ipot.lt.1 .or. ipot.gt.5) then
 740         write (iout,'(2a)') 'Error while reading SC interaction',
 741      &               'potential file - unknown potential type.'
 742 #ifdef MPI
 743         call MPI_Finalize(Ierror)
 744 #endif
 745         stop
 746       endif
 747       expon2=expon/2
 748       if(me.eq.king)
 749      & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
 750      & ', exponents are ',expon,2*expon
 751       goto (10,20,30,30,40) ipot
 752 C----------------------- LJ potential ---------------------------------
 753    10 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
 754      &   (sigma0(i),i=1,ntyp)
 755       if (lprint) then
 756         write (iout,'(/a/)') 'Parameters of the LJ potential:'
 757         write (iout,'(a/)') 'The epsilon array:'
 758         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 759         write (iout,'(/a)') 'One-body parameters:'
 760         write (iout,'(a,4x,a)') 'residue','sigma'
 761         write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
 762       endif
 763       goto 50
 764 C----------------------- LJK potential --------------------------------
 765    20 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
 766      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
 767       if (lprint) then
 768         write (iout,'(/a/)') 'Parameters of the LJK potential:'
 769         write (iout,'(a/)') 'The epsilon array:'
 770         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 771         write (iout,'(/a)') 'One-body parameters:'
 772         write (iout,'(a,4x,2a)') 'residue','   sigma  ','    r0    '
 773         write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
 774      &        i=1,ntyp)
 775       endif
 776       goto 50
 777 C---------------------- GB or BP potential -----------------------------
 778    30 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
 779      &  (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip(i),i=1,ntyp),
 780      &  (alp(i),i=1,ntyp)
 781 C For the GB potential convert sigma'**2 into chi'
 782       if (ipot.eq.4) then
 783         do i=1,ntyp
 784           chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
 785         enddo
 786       endif
 787       if (lprint) then
 788         write (iout,'(/a/)') 'Parameters of the BP potential:'
 789         write (iout,'(a/)') 'The epsilon array:'
 790         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 791         write (iout,'(/a)') 'One-body parameters:'
 792         write (iout,'(a,4x,4a)') 'residue','   sigma  ','s||/s_|_^2',
 793      &       '    chip  ','    alph  '
 794         write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
 795      &                     chip(i),alp(i),i=1,ntyp)
 796       endif
 797       goto 50
 798 C--------------------- GBV potential -----------------------------------
 799    40 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
 800      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
 801      &  (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
 802       if (lprint) then
 803         write (iout,'(/a/)') 'Parameters of the GBV potential:'
 804         write (iout,'(a/)') 'The epsilon array:'
 805         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 806         write (iout,'(/a)') 'One-body parameters:'
 807         write (iout,'(a,4x,5a)') 'residue','   sigma  ','    r0    ',
 808      &      's||/s_|_^2','    chip  ','    alph  '
 809         write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
 810      &           sigii(i),chip(i),alp(i),i=1,ntyp)
 811       endif
 812    50 continue
 813       close (isidep)
 814 C-----------------------------------------------------------------------
 815 C Calculate the "working" parameters of SC interactions.
 816       do i=2,ntyp
 817         do j=1,i-1
 818           eps(i,j)=eps(j,i)
 819         enddo
 820       enddo
 821       do i=1,ntyp
 822         do j=i,ntyp
 823           sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
 824           sigma(j,i)=sigma(i,j)
 825           rs0(i,j)=dwa16*sigma(i,j)
 826           rs0(j,i)=rs0(i,j)
 827         enddo
 828       enddo
 829       if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
 830      & 'Working parameters of the SC interactions:',
 831      & '     a    ','     b    ','   augm   ','  sigma ','   r0   ',
 832      & '  chi1   ','   chi2   '
 833       do i=1,ntyp
 834         do j=i,ntyp
 835           epsij=eps(i,j)
 836           if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
 837             rrij=sigma(i,j)
 838           else
 839             rrij=rr0(i)+rr0(j)
 840           endif
 841           r0(i,j)=rrij
 842           r0(j,i)=rrij
 843           rrij=rrij**expon
 844           epsij=eps(i,j)
 845           sigeps=dsign(1.0D0,epsij)
 846           epsij=dabs(epsij)
 847           aa(i,j)=epsij*rrij*rrij
 848           bb(i,j)=-sigeps*epsij*rrij
 849           aa(j,i)=aa(i,j)
 850           bb(j,i)=bb(i,j)
 851           if (ipot.gt.2) then
 852             sigt1sq=sigma0(i)**2
 853             sigt2sq=sigma0(j)**2
 854             sigii1=sigii(i)
 855             sigii2=sigii(j)
 856             ratsig1=sigt2sq/sigt1sq
 857             ratsig2=1.0D0/ratsig1
 858             chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
 859             if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
 860             rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
 861           else
 862             rsum_max=sigma(i,j)
 863           endif
 864 c         if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
 865             sigmaii(i,j)=rsum_max
 866             sigmaii(j,i)=rsum_max
 867 c         else
 868 c           sigmaii(i,j)=r0(i,j)
 869 c           sigmaii(j,i)=r0(i,j)
 870 c         endif
 871 cd        write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
 872           if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
 873             r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
 874             augm(i,j)=epsij*r_augm**(2*expon)
 875 c           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
 876             augm(j,i)=augm(i,j)
 877           else
 878             augm(i,j)=0.0D0
 879             augm(j,i)=0.0D0
 880           endif
 881           if (lprint) then
 882             write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
 883      &      restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
 884      &      sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
 885           endif
 886         enddo
 887       enddo
 888 #ifdef OLDSCP
 889 C
 890 C Define the SC-p interaction constants (hard-coded; old style)
 891 C
 892       do i=1,20
 893 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
 894 C helix formation)
 895 c       aad(i,1)=0.3D0*4.0D0**12
 896 C Following line for constants currently implemented
 897 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
 898         aad(i,1)=1.5D0*4.0D0**12
 899 c       aad(i,1)=0.17D0*5.6D0**12
 900         aad(i,2)=aad(i,1)
 901 C "Soft" SC-p repulsion
 902         bad(i,1)=0.0D0
 903 C Following line for constants currently implemented
 904 c       aad(i,1)=0.3D0*4.0D0**6
 905 C "Hard" SC-p repulsion
 906         bad(i,1)=3.0D0*4.0D0**6
 907 c       bad(i,1)=-2.0D0*0.17D0*5.6D0**6
 908         bad(i,2)=bad(i,1)
 909 c       aad(i,1)=0.0D0
 910 c       aad(i,2)=0.0D0
 911 c       bad(i,1)=1228.8D0
 912 c       bad(i,2)=1228.8D0
 913       enddo
 914 #else
 915 C
 916 C 8/9/01 Read the SC-p interaction constants from file
 917 C
 918       do i=1,ntyp
 919         read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
 920       enddo
 921       do i=1,ntyp
 922         aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
 923         aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
 924         bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
 925         bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
 926       enddo
 927
 928       if (lprint) then
 929         write (iout,*) "Parameters of SC-p interactions:"
 930         do i=1,20
 931           write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
 932      &     eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
 933         enddo
 934       endif
 935 #endif
 936 C
 937 C Define the constants of the disulfide bridge
 938 C
 939       ebr=-5.50D0
 940 c
 941 c Old arbitrary potential - commented out.
 942 c
 943 c      dbr= 4.20D0
 944 c      fbr= 3.30D0
 945 c
 946 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
 947 c energy surface of diethyl disulfide.
 948 c A. Liwo and U. Kozlowska, 11/24/03
 949 c
 950       D0CM = 3.78d0
 951       AKCM = 15.1d0
 952       AKTH = 11.0d0
 953       AKCT = 12.0d0
 954       V1SS =-1.08d0
 955       V2SS = 7.61d0
 956       V3SS = 13.7d0
 957 c      akcm=0.0d0
 958 c      akth=0.0d0
 959 c      akct=0.0d0
 960 c      v1ss=0.0d0
 961 c      v2ss=0.0d0
 962 c      v3ss=0.0d0
 963
 964       if(me.eq.king) then
 965       write (iout,'(/a)') "Disulfide bridge parameters:"
 966       write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
 967       write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
 968       write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
 969       write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
 970      &  ' v3ss:',v3ss
 971       endif
 972       return
 973   111 write (iout,*) "Error reading bending energy parameters."
 974       goto 999
 975   112 write (iout,*) "Error reading rotamer energy parameters."
 976       goto 999
 977   113 write (iout,*) "Error reading torsional energy parameters."
 978       goto 999
 979   114 write (iout,*) "Error reading double torsional energy parameters."
 980       goto 999
 981   115 write (iout,*)
 982      &  "Error reading cumulant (multibody energy) parameters."
 983       goto 999
 984   116 write (iout,*) "Error reading electrostatic energy parameters."
 985       goto 999
 986   117 write (iout,*) "Error reading side chain interaction parameters."
 987       goto 999
 988   118 write (iout,*) "Error reading SCp interaction parameters."
 989       goto 999
 990   119 write (iout,*) "Error reading SCCOR parameters"
 991   999 continue
 992 #ifdef MPI
 993       call MPI_Finalize(Ierror)
 994 #endif
 995       stop
 996       return
 997       end
 998
 999
1000       subroutine getenv_loc(var, val)
1001       character(*) var, val
1002
1003 #ifdef WINIFL
1004       character(2000) line
1005       external ilen
1006
1007       open (196,file='env',status='old',readonly,shared)
1008       iread=0
1009 c      write(*,*)'looking for ',var
1010 10    read(196,*,err=11,end=11)line
1011       iread=index(line,var)
1012 c      write(*,*)iread,' ',var,' ',line
1013       if (iread.eq.0) go to 10
1014 c      write(*,*)'---> ',line
1015 11    continue
1016       if(iread.eq.0) then
1017 c       write(*,*)'CHUJ'
1018        val=''
1019       else
1020        iread=iread+ilen(var)+1
1021        read (line(iread:),*,err=12,end=12) val
1022 c       write(*,*)'OK: ',var,' = ',val
1023       endif
1024       close(196)
1025       return
1026 12    val=''
1027       close(196)
1028 #elif (defined CRAY)
1029       integer lennam,lenval,ierror
1030 c
1031 c        getenv using a POSIX call, useful on the T3D
1032 c        Sept 1996, comment out error check on advice of H. Pritchard
1033 c
1034       lennam = len(var)
1035       if(lennam.le.0) stop '--error calling getenv--'
1036       call pxfgetenv(var,lennam,val,lenval,ierror)
1037 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1038 #else
1039       call getenv(var,val)
1040 #endif
1041
1042       return
1043       end