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