source/unres/src_MD/parmread.F

   1       subroutine parmread
   2 C
   3 C Read the parameters of the probability distributions of the virtual-bond
   4 C valence angles and the side chains and energy parameters.
   5 C
   6 C Important! Energy-term weights ARE NOT read here; they are read from the
   7 C main input file instead, because NO defaults have yet been set for these
   8 C parameters.
   9 C
  10       implicit real*8 (a-h,o-z)
  11       include 'DIMENSIONS'
  12 #ifdef MPI
  13       include "mpif.h"
  14       integer IERROR
  15 #endif
  16       include 'COMMON.IOUNITS'
  17       include 'COMMON.CHAIN'
  18       include 'COMMON.INTERACT'
  19       include 'COMMON.GEO'
  20       include 'COMMON.LOCAL'
  21       include 'COMMON.TORSION'
  22       include 'COMMON.SCCOR'
  23       include 'COMMON.SCROT'
  24       include 'COMMON.FFIELD'
  25       include 'COMMON.NAMES'
  26       include 'COMMON.SBRIDGE'
  27       include 'COMMON.MD'
  28       include 'COMMON.SETUP'
  29       character*1 t1,t2,t3
  30       character*1 onelett(4) /"G","A","P","D"/
  31       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 C      print *, "tu"
 214       nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
 215       write(iout,*) "I am here",ntyp1
 216       read (ithep,*,err=111,end=111) (ithetyp(i),i=1,ntyp1)
 217 C      write(iout,*) "I am herew"
 218       do i=-ntyp1,-1
 219         ithetyp(i)=-ithetyp(-i)
 220       enddo
 221       do iblock=1,2
 222       do i=-maxthetyp,maxthetyp
 223         do j=-maxthetyp,maxthetyp
 224           do k=-maxthetyp,maxthetyp
 225             aa0thet(i,j,k,iblock)=0.0d0
 226             do l=1,ntheterm
 227               aathet(l,i,j,k,iblock)=0.0d0
 228             enddo
 229             do l=1,ntheterm2
 230               do m=1,nsingle
 231                 bbthet(m,l,i,j,k,iblock)=0.0d0
 232                 ccthet(m,l,i,j,k,iblock)=0.0d0
 233                 ddthet(m,l,i,j,k,iblock)=0.0d0
 234                 eethet(m,l,i,j,k,iblock)=0.0d0
 235               enddo
 236             enddo
 237             do l=1,ntheterm3
 238               do m=1,ndouble
 239                 do mm=1,ndouble
 240                  ffthet(mm,m,l,i,j,k,iblock)=0.0d0
 241                  ggthet(mm,m,l,i,j,k,iblock)=0.0d0
 242                 enddo
 243               enddo
 244             enddo
 245           enddo
 246         enddo
 247       enddo
 248       enddo
 249 c VAR:iblock means terminally blocking group 1=non-proline 2=proline
 250       do iblock=1,2
 251 c VAR:ntethtyp is type of theta potentials type currently 0=glycine
 252 c VAR:1=non-glicyne non-proline 2=proline
 253 c VAR:negative values for D-aminoacid
 254       do i=0,nthetyp
 255 C        print *,i
 256         do j=-nthetyp,nthetyp
 257           do k=-nthetyp,nthetyp
 258             read (ithep,'(6a)',end=111,err=111) res1
 259 c VAR: aa0thet is variable describing the average value of Foureir
 260 c VAR: expansion series
 261             read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock)
 262 c VAR: aathet is foureir expansion in theta/2 angle for full formula
 263 c VAR: look at the fitting equation in Kozlowska et al., J. Phys.:
 264 Condens. Matter 19 (2007) 285203 and Sieradzan et al., unpublished
 265             read (ithep,*,end=111,err=111)
 266      &(aathet(l,i,j,k,iblock),l=1,ntheterm)
 267             read (ithep,*,end=111,err=111)
 268      &       ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
 269      &        (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
 270      &        (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
 271      &        (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle),
 272      &        ll=1,ntheterm2)
 273             read (ithep,*,end=111,err=111)
 274      &      (((ffthet(llll,lll,ll,i,j,k,iblock),
 275      &         ffthet(lll,llll,ll,i,j,k,iblock),
 276      &         ggthet(llll,lll,ll,i,j,k,iblock),
 277      &         ggthet(lll,llll,ll,i,j,k,iblock),
 278      &         llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
 279           enddo
 280         enddo
 281       enddo
 282
 283
 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 i=1,nthetyp+1
 362           do j=1,nthetyp+1
 363             do k=1,nthetyp+1
 364               write (iout,'(//4a)')
 365      &         'Type ',onelett(i),onelett(j),onelett(k)
 366               write (iout,'(//a,10x,a)') " l","a[l]"
 367               write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
 368               write (iout,'(i2,1pe15.5)')
 369      &           (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
 370             do l=1,ntheterm2
 371               write (iout,'(//2h m,4(9x,a,3h[m,,i1,1h]))')
 372      &          "b",l,"c",l,"d",l,"e",l
 373               do m=1,nsingle
 374                 write (iout,'(i2,4(1pe15.5))') m,
 375      &          bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),
 376      &          ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
 377               enddo
 378             enddo
 379             do l=1,ntheterm3
 380               write (iout,'(//3hm,n,4(6x,a,5h[m,n,,i1,1h]))')
 381      &          "f+",l,"f-",l,"g+",l,"g-",l
 382               do m=2,ndouble
 383                 do n=1,m-1
 384                   write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
 385      &              ffthet(n,m,l,i,j,k,iblock),
 386      &              ffthet(m,n,l,i,j,k,iblock),
 387      &              ggthet(n,m,l,i,j,k,iblock),
 388      &              ggthet(m,n,l,i,j,k,iblock)
 389                 enddo
 390               enddo
 391             enddo
 392           enddo
 393         enddo
 394       enddo
 395       call flush(iout)
 396       endif
 397       write (2,*) "Start reading THETA_PDB"
 398       do i=1,ntyp
 399         read (ithep_pdb,*,err=111,end=111)
 400      &     a0thet(i),(athet(j,i,1,1),j=1,2),
 401      &    (bthet(j,i,1,1),j=1,2)
 402         read (ithep_pdb,*,err=111,end=111) (polthet(j,i),j=0,3)
 403         read (ithep_pdb,*,err=111,end=111) (gthet(j,i),j=1,3)
 404         read (ithep_pdb,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
 405         sigc0(i)=sigc0(i)**2
 406       enddo
 407       do i=1,ntyp
 408       athet(1,i,1,-1)=athet(1,i,1,1)
 409       athet(2,i,1,-1)=athet(2,i,1,1)
 410       bthet(1,i,1,-1)=-bthet(1,i,1,1)
 411       bthet(2,i,1,-1)=-bthet(2,i,1,1)
 412       athet(1,i,-1,1)=-athet(1,i,1,1)
 413       athet(2,i,-1,1)=-athet(2,i,1,1)
 414       bthet(1,i,-1,1)=bthet(1,i,1,1)
 415       bthet(2,i,-1,1)=bthet(2,i,1,1)
 416       enddo
 417       do i=-ntyp,-1
 418       a0thet(i)=a0thet(-i)
 419       athet(1,i,-1,-1)=athet(1,-i,1,1)
 420       athet(2,i,-1,-1)=-athet(2,-i,1,1)
 421       bthet(1,i,-1,-1)=bthet(1,-i,1,1)
 422       bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
 423       athet(1,i,-1,1)=athet(1,-i,1,1)
 424       athet(2,i,-1,1)=-athet(2,-i,1,1)
 425       bthet(1,i,-1,1)=-bthet(1,-i,1,1)
 426       bthet(2,i,-1,1)=bthet(2,-i,1,1)
 427       athet(1,i,1,-1)=-athet(1,-i,1,1)
 428       athet(2,i,1,-1)=athet(2,-i,1,1)
 429       bthet(1,i,1,-1)=bthet(1,-i,1,1)
 430       bthet(2,i,1,-1)=-bthet(2,-i,1,1)
 431       theta0(i)=theta0(-i)
 432       sig0(i)=sig0(-i)
 433       sigc0(i)=sigc0(-i)
 434        do j=0,3
 435         polthet(j,i)=polthet(j,-i)
 436        enddo
 437        do j=1,3
 438          gthet(j,i)=gthet(j,-i)
 439        enddo
 440       enddo
 441       write (2,*) "End reading THETA_PDB"
 442       close (ithep_pdb)
 443 #endif
 444       close(ithep)
 445 #ifdef CRYST_SC
 446 C
 447 C Read the parameters of the probability distribution/energy expression
 448 C of the side chains.
 449 C
 450       do i=1,ntyp
 451         read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
 452         if (i.eq.10) then
 453           dsc_inv(i)=0.0D0
 454         else
 455           dsc_inv(i)=1.0D0/dsc(i)
 456         endif
 457         if (i.ne.10) then
 458         do j=1,nlob(i)
 459           do k=1,3
 460             do l=1,3
 461               blower(l,k,j)=0.0D0
 462             enddo
 463           enddo
 464         enddo
 465         bsc(1,i)=0.0D0
 466         read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),
 467      &    ((blower(k,l,1),l=1,k),k=1,3)
 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         enddo
 473         do j=1,nlob(i)
 474           do k=1,3
 475             do l=1,k
 476               akl=0.0D0
 477               do m=1,3
 478                 akl=akl+blower(k,m,j)*blower(l,m,j)
 479               enddo
 480               gaussc(k,l,j,i)=akl
 481               gaussc(l,k,j,i)=akl
 482             enddo
 483           enddo
 484         enddo
 485         endif
 486       enddo
 487       close (irotam)
 488       if (lprint) then
 489         write (iout,'(/a)') 'Parameters of side-chain local geometry'
 490         do i=1,ntyp
 491           nlobi=nlob(i)
 492           if (nlobi.gt.0) then
 493             if (LaTeX) then
 494               write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
 495      &         ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
 496                write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
 497      &                             'log h',(bsc(j,i),j=1,nlobi)
 498                write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
 499      &        'x',((censc(k,j,i),k=1,3),j=1,nlobi)
 500               do k=1,3
 501                 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
 502      &                 ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
 503               enddo
 504             else
 505               write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
 506               write (iout,'(a,f10.4,4(16x,f10.4))')
 507      &                             'Center  ',(bsc(j,i),j=1,nlobi)
 508               write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),
 509      &           j=1,nlobi)
 510               write (iout,'(a)')
 511             endif
 512           endif
 513         enddo
 514       endif
 515 #else
 516 C
 517 C Read scrot parameters for potentials determined from all-atom AM1 calculations
 518 C added by Urszula Kozlowska 07/11/2007
 519 C
 520       do i=1,ntyp
 521         read (irotam,*,end=112,err=112)
 522        if (i.eq.10) then
 523          read (irotam,*,end=112,err=112)
 524        else
 525          do j=1,65
 526            read(irotam,*,end=112,err=112) sc_parmin(j,i)
 527          enddo
 528        endif
 529       enddo
 530 C
 531 C Read the parameters of the probability distribution/energy expression
 532 C of the side chains.
 533 C
 534       write (2,*) "Start reading ROTAM_PDB"
 535
 536       do i=1,ntyp
 537         read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
 538         if (i.eq.10) then
 539           dsc_inv(i)=0.0D0
 540         else
 541           dsc_inv(i)=1.0D0/dsc(i)
 542         endif
 543         if (i.ne.10) then
 544         do j=1,nlob(i)
 545           do k=1,3
 546             do l=1,3
 547               blower(l,k,j)=0.0D0
 548             enddo
 549           enddo
 550         enddo
 551         bsc(1,i)=0.0D0
 552         read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),
 553      &    ((blower(k,l,1),l=1,k),k=1,3)
 554         do j=2,nlob(i)
 555           read (irotam_pdb,*,end=112,err=112) bsc(j,i)
 556           read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),
 557      &                                 ((blower(k,l,j),l=1,k),k=1,3)
 558         enddo
 559         do j=1,nlob(i)
 560           do k=1,3
 561             do l=1,k
 562               akl=0.0D0
 563               do m=1,3
 564                 akl=akl+blower(k,m,j)*blower(l,m,j)
 565               enddo
 566               gaussc(k,l,j,i)=akl
 567               gaussc(l,k,j,i)=akl
 568             enddo
 569           enddo
 570         enddo
 571         endif
 572       enddo
 573       close (irotam_pdb)
 574       write (2,*) "Ending reading ROTAM_PDB"
 575 #endif
 576       close(irotam)
 577
 578 #ifdef CRYST_TOR
 579 C
 580 C Read torsional parameters in old format
 581 C
 582       read (itorp,*,end=113,err=113) ntortyp,nterm_old
 583       if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
 584       read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
 585       do i=1,ntortyp
 586         do j=1,ntortyp
 587           read (itorp,'(a)')
 588           do k=1,nterm_old
 589             read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
 590           enddo
 591         enddo
 592       enddo
 593       close (itorp)
 594       if (lprint) then
 595         write (iout,'(/a/)') 'Torsional constants:'
 596         do i=1,ntortyp
 597           do j=1,ntortyp
 598             write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
 599             write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
 600           enddo
 601         enddo
 602       endif
 603 #else
 604 C
 605 C Read torsional parameters
 606 C
 607       read (itorp,*,end=113,err=113) ntortyp
 608       read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
 609       do iblock=1,2
 610       do i=-ntyp,-1
 611        itortyp(i)=-itortyp(-i)
 612       enddo
 613       write (iout,*) 'ntortyp',ntortyp
 614       do i=0,ntortyp-1
 615         do j=-ntortyp+1,ntortyp-1
 616           read (itorp,*,end=113,err=113) nterm(i,j,iblock),
 617      &          nlor(i,j,iblock)
 618           nterm(-i,-j,iblock)=nterm(i,j,iblock)
 619           nlor(-i,-j,iblock)=nlor(i,j,iblock)
 620           v0ij=0.0d0
 621           si=-1.0d0
 622           do k=1,nterm(i,j,iblock)
 623             read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),
 624      &      v2(k,i,j,iblock)
 625             v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
 626             v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
 627             v0ij=v0ij+si*v1(k,i,j,iblock)
 628             si=-si
 629 c         write(iout,*) i,j,k,iblock,nterm(i,j,iblock)
 630 c         write(iout,*) v1(k,-i,-j,iblock),v1(k,i,j,iblock),
 631 c      &v2(k,-i,-j,iblock),v2(k,i,j,iblock)
 632           enddo
 633           do k=1,nlor(i,j,iblock)
 634             read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),
 635      &        vlor2(k,i,j),vlor3(k,i,j)
 636             v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
 637           enddo
 638           v0(i,j,iblock)=v0ij
 639           v0(-i,-j,iblock)=v0ij
 640         enddo
 641       enddo
 642       enddo
 643       close (itorp)
 644       if (lprint) then
 645         write (iout,'(/a/)') 'Torsional constants:'
 646         do i=1,ntortyp
 647           do j=1,ntortyp
 648             write (iout,*) 'ityp',i,' jtyp',j
 649             write (iout,*) 'Fourier constants'
 650             do k=1,nterm(i,j,iblock)
 651               write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
 652      &        v2(k,i,j,iblock)
 653             enddo
 654             write (iout,*) 'Lorenz constants'
 655             do k=1,nlor(i,j,iblock)
 656               write (iout,'(3(1pe15.5))')
 657      &         vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
 658             enddo
 659           enddo
 660         enddo
 661       endif
 662 C
 663 C 6/23/01 Read parameters for double torsionals
 664 C
 665 C      do i=1,ntortyp
 666 C        do j=1,ntortyp
 667 C          do k=1,ntortyp
 668       do iblock=1,2
 669       do i=0,ntortyp-1
 670         do j=-ntortyp+1,ntortyp-1
 671           do k=-ntortyp+1,ntortyp-1
 672             read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
 673               write (iout,*) "OK onelett",
 674      &         i,j,k,t1,t2,t3
 675
 676             if (t1.ne.onelett(i) .or. t2.ne.onelett(j)
 677      &        .or. t3.ne.onelett(k)) then
 678               write (iout,*) "Error in double torsional parameter file",
 679      &         i,j,k,t1,t2,t3
 680 #ifdef MPI
 681               call MPI_Finalize(Ierror)
 682 #endif
 683                stop "Error in double torsional parameter file"
 684             endif
 685            read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),
 686      &         ntermd_2(i,j,k,iblock)
 687             ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
 688             ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
 689             read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,
 690      &         ntermd_1(i,j,k,iblock))
 691             read (itordp,*,end=114,err=114)(v1s(1,l,i,j,k,iblock),l=1,
 692      &         ntermd_1(i,j,k,iblock))
 693             read (itordp,*,end=114,err=114)(v1c(2,l,i,j,k,iblock),l=1,
 694      &         ntermd_1(i,j,k,iblock))
 695             read (itordp,*,end=114,err=114)(v1s(2,l,i,j,k,iblock),l=1,
 696      &         ntermd_1(i,j,k,iblock))
 697 C Martix of D parameters for one dimesional foureir series
 698             do l=1,ntermd_1(i,j,k,iblock)
 699              v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
 700              v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
 701              v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
 702              v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
 703 c            write(iout,*) "whcodze" ,
 704 c     & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
 705             enddo
 706             read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),
 707      &         v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
 708      &         v2s(m,l,i,j,k,iblock),
 709      &         m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
 710 C Martix of D parameters for two dimesional fourier series
 711             do l=1,ntermd_2(i,j,k,iblock)
 712              do m=1,l-1
 713              v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
 714              v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
 715              v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
 716              v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
 717              enddo!m
 718             enddo!l
 719           enddo!k
 720         enddo!j
 721       enddo!i
 722       enddo!iblock
 723       if (lprint) then
 724       write (iout,*)
 725       write (iout,*) 'Constants for double torsionals'
 726       do iblock=1,2
 727       do i=0,ntortyp-1
 728         do j=-ntortyp+1,ntortyp-1
 729           do k=-ntortyp+1,ntortyp-1
 730             write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
 731      &        ' nsingle',ntermd_1(i,j,k,iblock),
 732      &        ' ndouble',ntermd_2(i,j,k,iblock)
 733             write (iout,*)
 734             write (iout,*) 'Single angles:'
 735             do l=1,ntermd_1(i,j,k,iblock)
 736               write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
 737      &           v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
 738      &           v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),
 739      &           v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
 740             enddo
 741             write (iout,*)
 742             write (iout,*) 'Pairs of angles:'
 743             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
 744             do l=1,ntermd_2(i,j,k,iblock)
 745               write (iout,'(i5,20f10.5)')
 746      &         l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
 747             enddo
 748             write (iout,*)
 749             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
 750             do l=1,ntermd_2(i,j,k,iblock)
 751               write (iout,'(i5,20f10.5)')
 752      &         l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),
 753      &         (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
 754             enddo
 755             write (iout,*)
 756           enddo
 757         enddo
 758       enddo
 759       enddo
 760       endif
 761 #endif
 762 C Read of Side-chain backbone correlation parameters
 763 C Modified 11 May 2012 by Adasko
 764 CCC
 765 C
 766       read (isccor,*,end=1113,err=1113) nsccortyp
 767 #ifdef SCCORPDB
 768       write (iout,*) "Tu wchodze"
 769       read (isccor,*,end=1113,err=1113) (isccortyp(i),i=1,ntyp)
 770       do i=-ntyp,-1
 771         isccortyp(i)=-isccortyp(-i)
 772       enddo
 773       iscprol=isccortyp(20)
 774 C      write (iout,*) 'ntortyp',ntortyp
 775       maxinter=3
 776 cc maxinter is maximum interaction sites
 777       do l=1,maxinter
 778       do i=1,nsccortyp
 779         do j=1,nsccortyp
 780          read (isccor,*,end=1113,err=1113) nterm_sccor(i,j),
 781      &             nlor_sccor(i,j)
 782           print *,i,j,l
 783           v0ijsccor=0.0d0
 784           v0ijsccor1=0.0d0
 785           v0ijsccor2=0.0d0
 786           v0ijsccor3=0.0d0
 787           si=-1.0d0
 788           nterm_sccor(-i,j)=nterm_sccor(i,j)
 789           nterm_sccor(-i,-j)=nterm_sccor(i,j)
 790           nterm_sccor(i,-j)=nterm_sccor(i,j)
 791           do k=1,nterm_sccor(i,j)
 792            read (isccor,*,end=1113,err=1113) kk,v1sccor(k,l,i,j)
 793      &    ,v2sccor(k,l,i,j)
 794 c            write(iout,*) "k=",kk
 795             if (j.eq.iscprol) then
 796               if (i.eq.isccortyp(10)) then
 797               v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
 798               v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
 799               else
 800              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
 801      &                        +v2sccor(k,l,i,j)*dsqrt(0.75d0)
 802              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
 803      &                        +v1sccor(k,l,i,j)*dsqrt(0.75d0)
 804              v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
 805              v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
 806              v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
 807              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
 808              endif
 809             else
 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                 if (j.eq.isccortyp(10)) then
 815               v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
 816               v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
 817                 else
 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              v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
 823              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
 824             endif
 825              endif
 826              endif
 827 C         read (isccor,*,end=1113,err=1113) kk,v1sccor(k,l,i,j)
 828 C     &    ,v2sccor(k,l,i,j)
 829             v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
 830             v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
 831             v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
 832             v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
 833             si=-si
 834           enddo
 835          do k=1,nlor_sccor(i,j)
 836             read (isccor,*,end=1113,err=1113) kk,vlor1sccor(k,i,j),
 837      &        vlor2sccor(k,i,j),vlor3sccor(k,i,j)
 838             v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
 839      &(1+vlor3sccor(k,i,j)**2)
 840           enddo
 841           v0sccor(l,i,j)=v0ijsccor
 842           v0sccor(l,-i,j)=v0ijsccor1
 843           v0sccor(l,i,-j)=v0ijsccor2
 844           v0sccor(l,-i,-j)=v0ijsccor3
 845         enddo
 846       enddo
 847       enddo
 848       close (isccor)
 849 #else
 850       write(iout,*) "a tu nie wchodze"
 851       read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
 852 c      write (iout,*) 'ntortyp',ntortyp
 853       maxinter=3
 854 cc maxinter is maximum interaction sites
 855       do l=1,maxinter
 856       do i=1,nsccortyp
 857         do j=1,nsccortyp
 858           read (isccor,*,end=119,err=119)
 859      & nterm_sccor(i,j),nlor_sccor(i,j)
 860           v0ijsccor=0.0d0
 861           si=-1.0d0
 862
 863           do k=1,nterm_sccor(i,j)
 864             read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
 865      &    ,v2sccor(k,l,i,j)
 866             v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
 867             si=-si
 868           enddo
 869           do k=1,nlor_sccor(i,j)
 870             read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
 871      &        vlor2sccor(k,i,j),vlor3sccor(k,i,j)
 872             v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
 873      &(1+vlor3sccor(k,i,j)**2)
 874           enddo
 875           v0sccor(i,j,iblock)=v0ijsccor
 876         enddo
 877       enddo
 878       enddo
 879       close (isccor)
 880
 881 #endif
 882       if (lprint) then
 883         write (iout,'(/a/)') 'Torsional constants:'
 884         do i=1,nsccortyp
 885           do j=1,nsccortyp
 886             write (iout,*) 'ityp',i,' jtyp',j
 887             write (iout,*) 'Fourier constants'
 888             do k=1,nterm_sccor(i,j)
 889       write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
 890             enddo
 891             write (iout,*) 'Lorenz constants'
 892             do k=1,nlor_sccor(i,j)
 893               write (iout,'(3(1pe15.5))')
 894      &         vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
 895             enddo
 896           enddo
 897         enddo
 898       endif
 899 C
 900 C
 901 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
 902 C         interaction energy of the Gly, Ala, and Pro prototypes.
 903 C
 904       if (lprint) then
 905         write (iout,*)
 906         write (iout,*) "Coefficients of the cumulants"
 907       endif
 908       read (ifourier,*) nloctyp
 909       do i=1,nloctyp
 910         read (ifourier,*,end=115,err=115)
 911         read (ifourier,*,end=115,err=115) (b(ii),ii=1,13)
 912         if (lprint) then
 913         write (iout,*) 'Type',i
 914         write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
 915         endif
 916         B1(1,i)  = b(3)
 917         B1(2,i)  = b(5)
 918 c        b1(1,i)=0.0d0
 919 c        b1(2,i)=0.0d0
 920         B1tilde(1,i) = b(3)
 921         B1tilde(2,i) =-b(5)
 922 C        B1tilde(1,-i) =-b(3)
 923 C        B1tilde(2,-i) =b(5)
 924 c        b1tilde(1,i)=0.0d0
 925 c        b1tilde(2,i)=0.0d0
 926         B2(1,i)  = b(2)
 927         B2(2,i)  = b(4)
 928 c        b2(1,i)=0.0d0
 929 c        b2(2,i)=0.0d0
 930         CC(1,1,i)= b(7)
 931         CC(2,2,i)=-b(7)
 932         CC(2,1,i)= b(9)
 933         CC(1,2,i)= b(9)
 934 c        CC(1,1,i)=0.0d0
 935 c        CC(2,2,i)=0.0d0
 936 c        CC(2,1,i)=0.0d0
 937 c        CC(1,2,i)=0.0d0
 938         Ctilde(1,1,i)=b(7)
 939         Ctilde(1,2,i)=b(9)
 940         Ctilde(2,1,i)=-b(9)
 941         Ctilde(2,2,i)=b(7)
 942 c        Ctilde(1,1,i)=0.0d0
 943 c        Ctilde(1,2,i)=0.0d0
 944 c        Ctilde(2,1,i)=0.0d0
 945 c        Ctilde(2,2,i)=0.0d0
 946         DD(1,1,i)= b(6)
 947         DD(2,2,i)=-b(6)
 948         DD(2,1,i)= b(8)
 949         DD(1,2,i)= b(8)
 950 c        DD(1,1,i)=0.0d0
 951 c        DD(2,2,i)=0.0d0
 952 c        DD(2,1,i)=0.0d0
 953 c        DD(1,2,i)=0.0d0
 954         Dtilde(1,1,i)=b(6)
 955         Dtilde(1,2,i)=b(8)
 956         Dtilde(2,1,i)=-b(8)
 957         Dtilde(2,2,i)=b(6)
 958 c        Dtilde(1,1,i)=0.0d0
 959 c        Dtilde(1,2,i)=0.0d0
 960 c        Dtilde(2,1,i)=0.0d0
 961 c        Dtilde(2,2,i)=0.0d0
 962         EE(1,1,i)= b(10)+b(11)
 963         EE(2,2,i)=-b(10)+b(11)
 964         EE(2,1,i)= b(12)-b(13)
 965         EE(1,2,i)= b(12)+b(13)
 966 c        ee(1,1,i)=1.0d0
 967 c        ee(2,2,i)=1.0d0
 968 c        ee(2,1,i)=0.0d0
 969 c        ee(1,2,i)=0.0d0
 970 c        ee(2,1,i)=ee(1,2,i)
 971       enddo
 972       if (lprint) then
 973       do i=1,nloctyp
 974         write (iout,*) 'Type',i
 975         write (iout,*) 'B1'
 976         write(iout,*) B1(1,i),B1(2,i)
 977         write (iout,*) 'B2'
 978         write(iout,*) B2(1,i),B2(2,i)
 979         write (iout,*) 'CC'
 980         do j=1,2
 981           write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
 982         enddo
 983         write(iout,*) 'DD'
 984         do j=1,2
 985           write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
 986         enddo
 987         write(iout,*) 'EE'
 988         do j=1,2
 989           write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
 990         enddo
 991       enddo
 992       endif
 993 C
 994 C Read electrostatic-interaction parameters
 995 C
 996       if (lprint) then
 997         write (iout,*)
 998         write (iout,'(/a)') 'Electrostatic interaction constants:'
 999         write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
1000      &            'IT','JT','APP','BPP','AEL6','AEL3'
1001       endif
1002       read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
1003       read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
1004       read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
1005       read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
1006       close (ielep)
1007       do i=1,2
1008         do j=1,2
1009         rri=rpp(i,j)**6
1010         app (i,j)=epp(i,j)*rri*rri
1011         bpp (i,j)=-2.0D0*epp(i,j)*rri
1012         ael6(i,j)=elpp6(i,j)*4.2D0**6
1013         ael3(i,j)=elpp3(i,j)*4.2D0**3
1014         if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
1015      &                    ael6(i,j),ael3(i,j)
1016         enddo
1017       enddo
1018 C
1019 C Read side-chain interaction parameters.
1020 C
1021       read (isidep,*,end=117,err=117) ipot,expon
1022       if (ipot.lt.1 .or. ipot.gt.5) then
1023         write (iout,'(2a)') 'Error while reading SC interaction',
1024      &               'potential file - unknown potential type.'
1025 #ifdef MPI
1026         call MPI_Finalize(Ierror)
1027 #endif
1028         stop
1029       endif
1030       expon2=expon/2
1031       if(me.eq.king)
1032      & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
1033      & ', exponents are ',expon,2*expon
1034       goto (10,20,30,30,40) ipot
1035 C----------------------- LJ potential ---------------------------------
1036    10 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
1037      &   (sigma0(i),i=1,ntyp)
1038       if (lprint) then
1039         write (iout,'(/a/)') 'Parameters of the LJ potential:'
1040         write (iout,'(a/)') 'The epsilon array:'
1041         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1042         write (iout,'(/a)') 'One-body parameters:'
1043         write (iout,'(a,4x,a)') 'residue','sigma'
1044         write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1045       endif
1046       goto 50
1047 C----------------------- LJK potential --------------------------------
1048    20 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
1049      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1050       if (lprint) then
1051         write (iout,'(/a/)') 'Parameters of the LJK potential:'
1052         write (iout,'(a/)') 'The epsilon array:'
1053         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1054         write (iout,'(/a)') 'One-body parameters:'
1055         write (iout,'(a,4x,2a)') 'residue','   sigma  ','    r0    '
1056         write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
1057      &        i=1,ntyp)
1058       endif
1059       goto 50
1060 C---------------------- GB or BP potential -----------------------------
1061    30 do i=1,ntyp
1062         read (isidep,*,end=116,err=116)(eps(i,j),j=i,ntyp)
1063        enddo
1064        read (isidep,*,end=116,err=116)(sigma0(i),i=1,ntyp)
1065        read (isidep,*,end=116,err=116)(sigii(i),i=1,ntyp)
1066        read (isidep,*,end=116,err=116)(chip(i),i=1,ntyp)
1067        read (isidep,*,end=116,err=116)(alp(i),i=1,ntyp)
1068
1069 c   30 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
1070 c     &  (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip(i),i=1,ntyp),
1071 c     &  (alp(i),i=1,ntyp)
1072 C For the GB potential convert sigma'**2 into chi'
1073       if (ipot.eq.4) then
1074         do i=1,ntyp
1075           chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1076         enddo
1077       endif
1078       if (lprint) then
1079         write (iout,'(/a/)') 'Parameters of the BP potential:'
1080         write (iout,'(a/)') 'The epsilon array:'
1081         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1082         write (iout,'(/a)') 'One-body parameters:'
1083         write (iout,'(a,4x,4a)') 'residue','   sigma  ','s||/s_|_^2',
1084      &       '    chip  ','    alph  '
1085         write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1086      &                     chip(i),alp(i),i=1,ntyp)
1087       endif
1088       goto 50
1089 C--------------------- GBV potential -----------------------------------
1090    40 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
1091      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1092      &  (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1093       if (lprint) then
1094         write (iout,'(/a/)') 'Parameters of the GBV potential:'
1095         write (iout,'(a/)') 'The epsilon array:'
1096         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1097         write (iout,'(/a)') 'One-body parameters:'
1098         write (iout,'(a,4x,5a)') 'residue','   sigma  ','    r0    ',
1099      &      's||/s_|_^2','    chip  ','    alph  '
1100         write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1101      &           sigii(i),chip(i),alp(i),i=1,ntyp)
1102       endif
1103    50 continue
1104       close (isidep)
1105 C-----------------------------------------------------------------------
1106 C Calculate the "working" parameters of SC interactions.
1107       do i=2,ntyp
1108         do j=1,i-1
1109           eps(i,j)=eps(j,i)
1110         enddo
1111       enddo
1112       do i=1,ntyp
1113         do j=i,ntyp
1114           sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1115           sigma(j,i)=sigma(i,j)
1116           rs0(i,j)=dwa16*sigma(i,j)
1117           rs0(j,i)=rs0(i,j)
1118         enddo
1119       enddo
1120       if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1121      & 'Working parameters of the SC interactions:',
1122      & '     a    ','     b    ','   augm   ','  sigma ','   r0   ',
1123      & '  chi1   ','   chi2   '
1124       do i=1,ntyp
1125         do j=i,ntyp
1126           epsij=eps(i,j)
1127           if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1128             rrij=sigma(i,j)
1129           else
1130             rrij=rr0(i)+rr0(j)
1131           endif
1132           r0(i,j)=rrij
1133           r0(j,i)=rrij
1134           rrij=rrij**expon
1135           epsij=eps(i,j)
1136           sigeps=dsign(1.0D0,epsij)
1137           epsij=dabs(epsij)
1138           aa(i,j)=epsij*rrij*rrij
1139           bb(i,j)=-sigeps*epsij*rrij
1140           aa(j,i)=aa(i,j)
1141           bb(j,i)=bb(i,j)
1142           if (ipot.gt.2) then
1143             sigt1sq=sigma0(i)**2
1144             sigt2sq=sigma0(j)**2
1145             sigii1=sigii(i)
1146             sigii2=sigii(j)
1147             ratsig1=sigt2sq/sigt1sq
1148             ratsig2=1.0D0/ratsig1
1149             chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1150             if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1151             rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1152           else
1153             rsum_max=sigma(i,j)
1154           endif
1155 c         if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1156             sigmaii(i,j)=rsum_max
1157             sigmaii(j,i)=rsum_max
1158 c         else
1159 c           sigmaii(i,j)=r0(i,j)
1160 c           sigmaii(j,i)=r0(i,j)
1161 c         endif
1162 cd        write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1163           if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1164             r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1165             augm(i,j)=epsij*r_augm**(2*expon)
1166 c           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1167             augm(j,i)=augm(i,j)
1168           else
1169             augm(i,j)=0.0D0
1170             augm(j,i)=0.0D0
1171           endif
1172           if (lprint) then
1173             write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1174      &      restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1175      &      sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1176           endif
1177         enddo
1178       enddo
1179 #ifdef OLDSCP
1180 C
1181 C Define the SC-p interaction constants (hard-coded; old style)
1182 C
1183       do i=1,ntyp
1184 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1185 C helix formation)
1186 c       aad(i,1)=0.3D0*4.0D0**12
1187 C Following line for constants currently implemented
1188 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1189         aad(i,1)=1.5D0*4.0D0**12
1190 c       aad(i,1)=0.17D0*5.6D0**12
1191         aad(i,2)=aad(i,1)
1192 C "Soft" SC-p repulsion
1193         bad(i,1)=0.0D0
1194 C Following line for constants currently implemented
1195 c       aad(i,1)=0.3D0*4.0D0**6
1196 C "Hard" SC-p repulsion
1197         bad(i,1)=3.0D0*4.0D0**6
1198 c       bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1199         bad(i,2)=bad(i,1)
1200 c       aad(i,1)=0.0D0
1201 c       aad(i,2)=0.0D0
1202 c       bad(i,1)=1228.8D0
1203 c       bad(i,2)=1228.8D0
1204       enddo
1205 #else
1206 C
1207 C 8/9/01 Read the SC-p interaction constants from file
1208 C
1209       do i=1,ntyp
1210         read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1211       enddo
1212       do i=1,ntyp
1213         aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1214         aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1215         bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1216         bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1217       enddo
1218
1219       if (lprint) then
1220         write (iout,*) "Parameters of SC-p interactions:"
1221         do i=1,ntyp
1222           write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1223      &     eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1224         enddo
1225       endif
1226 #endif
1227 C
1228 C Define the constants of the disulfide bridge
1229 C
1230       ebr=-5.50D0
1231 c
1232 c Old arbitrary potential - commented out.
1233 c
1234 c      dbr= 4.20D0
1235 c      fbr= 3.30D0
1236 c
1237 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1238 c energy surface of diethyl disulfide.
1239 c A. Liwo and U. Kozlowska, 11/24/03
1240 c
1241       D0CM = 3.78d0
1242       AKCM = 15.1d0
1243       AKTH = 11.0d0
1244       AKCT = 12.0d0
1245       V1SS =-1.08d0
1246       V2SS = 7.61d0
1247       V3SS = 13.7d0
1248 c      akcm=0.0d0
1249 c      akth=0.0d0
1250 c      akct=0.0d0
1251 c      v1ss=0.0d0
1252 c      v2ss=0.0d0
1253 c      v3ss=0.0d0
1254
1255       if(me.eq.king) then
1256       write (iout,'(/a)') "Disulfide bridge parameters:"
1257       write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1258       write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1259       write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1260       write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1261      &  ' v3ss:',v3ss
1262       endif
1263       return
1264   111 write (iout,*) "Error reading bending energy parameters."
1265       goto 999
1266   112 write (iout,*) "Error reading rotamer energy parameters."
1267       goto 999
1268   113 write (iout,*) "Error reading torsional energy parameters."
1269       goto 999
1270  1113 write (iout,*)
1271      &  "Error reading side-chain torsional energy parameters."
1272       goto 999
1273   114 write (iout,*) "Error reading double torsional energy parameters."
1274       goto 999
1275   115 write (iout,*)
1276      &  "Error reading cumulant (multibody energy) parameters."
1277       goto 999
1278   116 write (iout,*) "Error reading electrostatic energy parameters."
1279       goto 999
1280   117 write (iout,*) "Error reading side chain interaction parameters."
1281       goto 999
1282   118 write (iout,*) "Error reading SCp interaction parameters."
1283       goto 999
1284   119 write (iout,*) "Error reading SCCOR parameters"
1285   999 continue
1286 #ifdef MPI
1287       call MPI_Finalize(Ierror)
1288 #endif
1289       stop
1290       return
1291       end
1292
1293
1294       subroutine getenv_loc(var, val)
1295       character(*) var, val
1296
1297 #ifdef WINIFL
1298       character(2000) line
1299       external ilen
1300
1301       open (196,file='env',status='old',readonly,shared)
1302       iread=0
1303 c      write(*,*)'looking for ',var
1304 10    read(196,*,err=11,end=11)line
1305       iread=index(line,var)
1306 c      write(*,*)iread,' ',var,' ',line
1307       if (iread.eq.0) go to 10
1308 c      write(*,*)'---> ',line
1309 11    continue
1310       if(iread.eq.0) then
1311 c       write(*,*)'CHUJ'
1312        val=''
1313       else
1314        iread=iread+ilen(var)+1
1315        read (line(iread:),*,err=12,end=12) val
1316 c       write(*,*)'OK: ',var,' = ',val
1317       endif
1318       close(196)
1319       return
1320 12    val=''
1321       close(196)
1322 #elif (defined CRAY)
1323       integer lennam,lenval,ierror
1324 c
1325 c        getenv using a POSIX call, useful on the T3D
1326 c        Sept 1996, comment out error check on advice of H. Pritchard
1327 c
1328       lennam = len(var)
1329       if(lennam.le.0) stop '--error calling getenv--'
1330       call pxfgetenv(var,lennam,val,lenval,ierror)
1331 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1332 #else
1333       call getenv(var,val)
1334 #endif
1335
1336       return
1337       end