source/unres/src_MD-M-newcorr/parmread.F

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