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