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