source/unres/src_MD-M/parmread.F

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