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 of Side-chain backbone correlation parameters
 779 C Modified 11 May 2012 by Adasko
 780 CCC
 781 C
 782       read (isccor,*,end=119,err=119) nsccortyp
 783 #ifdef SCCORPDB
 784       read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
 785       do i=-ntyp,-1
 786         isccortyp(i)=-isccortyp(-i)
 787       enddo
 788       iscprol=isccortyp(20)
 789 c      write (iout,*) 'ntortyp',ntortyp
 790       maxinter=3
 791 cc maxinter is maximum interaction sites
 792       do l=1,maxinter
 793       do i=1,nsccortyp
 794         do j=1,nsccortyp
 795           read (isccor,*,end=119,err=119)
 796      &nterm_sccor(i,j),nlor_sccor(i,j)
 797           v0ijsccor=0.0d0
 798           v0ijsccor1=0.0d0
 799           v0ijsccor2=0.0d0
 800           v0ijsccor3=0.0d0
 801           si=-1.0d0
 802           nterm_sccor(-i,j)=nterm_sccor(i,j)
 803           nterm_sccor(-i,-j)=nterm_sccor(i,j)
 804           nterm_sccor(i,-j)=nterm_sccor(i,j)
 805           do k=1,nterm_sccor(i,j)
 806             read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
 807      &    ,v2sccor(k,l,i,j)
 808             if (j.eq.iscprol) then
 809              if (i.eq.isccortyp(10)) then
 810              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
 811              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
 812              else
 813              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
 814      &                        +v2sccor(k,l,i,j)*dsqrt(0.75d0)
 815              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
 816      &                        +v1sccor(k,l,i,j)*dsqrt(0.75d0)
 817              v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
 818              v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
 819              v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
 820              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
 821              endif
 822             else
 823              if (i.eq.isccortyp(10)) then
 824              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
 825              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
 826              else
 827                if (j.eq.isccortyp(10)) then
 828              v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
 829              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
 830                else
 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              v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
 836              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
 837                 endif
 838                endif
 839             endif
 840             v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
 841             v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
 842             v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
 843             v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
 844             si=-si
 845           enddo
 846           do k=1,nlor_sccor(i,j)
 847             read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
 848      &        vlor2sccor(k,i,j),vlor3sccor(k,i,j)
 849             v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
 850      &(1+vlor3sccor(k,i,j)**2)
 851           enddo
 852           v0sccor(l,i,j)=v0ijsccor
 853           v0sccor(l,-i,j)=v0ijsccor1
 854           v0sccor(l,i,-j)=v0ijsccor2
 855           v0sccor(l,-i,-j)=v0ijsccor3
 856         enddo
 857       enddo
 858       enddo
 859       close (isccor)
 860 #else
 861       read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
 862 c      write (iout,*) 'ntortyp',ntortyp
 863       maxinter=3
 864 cc maxinter is maximum interaction sites
 865       do l=1,maxinter
 866       do i=1,nsccortyp
 867         do j=1,nsccortyp
 868           read (isccor,*,end=119,err=119)
 869      & nterm_sccor(i,j),nlor_sccor(i,j)
 870           v0ijsccor=0.0d0
 871           si=-1.0d0
 872
 873           do k=1,nterm_sccor(i,j)
 874             read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
 875      &    ,v2sccor(k,l,i,j)
 876             v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
 877             si=-si
 878           enddo
 879           do k=1,nlor_sccor(i,j)
 880             read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
 881      &        vlor2sccor(k,i,j),vlor3sccor(k,i,j)
 882             v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
 883      &(1+vlor3sccor(k,i,j)**2)
 884           enddo
 885           v0sccor(i,j,iblock)=v0ijsccor
 886         enddo
 887       enddo
 888       enddo
 889       close (isccor)
 890
 891 #endif
 892       if (lprint) then
 893         write (iout,'(/a/)') 'Torsional constants:'
 894         do i=1,nsccortyp
 895           do j=1,nsccortyp
 896             write (iout,*) 'ityp',i,' jtyp',j
 897             write (iout,*) 'Fourier constants'
 898             do k=1,nterm_sccor(i,j)
 899       write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
 900             enddo
 901             write (iout,*) 'Lorenz constants'
 902             do k=1,nlor_sccor(i,j)
 903               write (iout,'(3(1pe15.5))')
 904      &         vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
 905             enddo
 906           enddo
 907         enddo
 908       endif
 909
 910 C
 911 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
 912 C         interaction energy of the Gly, Ala, and Pro prototypes.
 913 C
 914       if (lprint) then
 915         write (iout,*)
 916         write (iout,*) "Coefficients of the cumulants"
 917       endif
 918       read (ifourier,*) nloctyp
 919
 920       do i=0,nloctyp-1
 921         read (ifourier,*,end=115,err=115)
 922         read (ifourier,*,end=115,err=115) (b(ii,i),ii=1,13)
 923 #ifdef NEWCORR
 924         read (ifourier,*,end=115,err=115) (bnew1(ii,1,i),ii=1,3)
 925         read (ifourier,*,end=115,err=115) (bnew2(ii,1,i),ii=1,3)
 926         read (ifourier,*,end=115,err=115) (bnew1(ii,2,i),ii=1,1)
 927         read (ifourier,*,end=115,err=115) (bnew2(ii,2,i),ii=1,1)
 928         read (ifourier,*,end=115,err=115) (eenew(ii,i),ii=1,1)
 929 #endif
 930         if (lprint) then
 931         write (iout,*) 'Type',i
 932         write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
 933         endif
 934 c        B1(1,i)  = b(3)
 935 c        B1(2,i)  = b(5)
 936 c        B1(1,-i) = b(3)
 937 c        B1(2,-i) = -b(5)
 938 c        b1(1,i)=0.0d0
 939 c        b1(2,i)=0.0d0
 940 c        B1tilde(1,i) = b(3)
 941 c        B1tilde(2,i) =-b(5)
 942 c        B1tilde(1,-i) =-b(3)
 943 c        B1tilde(2,-i) =b(5)
 944 c        b1tilde(1,i)=0.0d0
 945 c        b1tilde(2,i)=0.0d0
 946 c        B2(1,i)  = b(2)
 947 c        B2(2,i)  = b(4)
 948 c        B2(1,-i)  =b(2)
 949 c        B2(2,-i)  =-b(4)
 950         B1tilde(1,i) = b(3,i)
 951         B1tilde(2,i) =-b(5,i)
 952         B1tilde(1,-i) =-b(3,i)
 953         B1tilde(2,-i) =b(5,i)
 954         b1tilde(1,i)=0.0d0
 955         b1tilde(2,i)=0.0d0
 956         B2(1,i)  = b(2,i)
 957         B2(2,i)  = b(4,i)
 958         B2(1,-i)  =b(2,i)
 959         B2(2,-i)  =-b(4,i)
 960
 961 c        b2(1,i)=0.0d0
 962 c        b2(2,i)=0.0d0
 963         CC(1,1,i)= b(7,i)
 964         CC(2,2,i)=-b(7,i)
 965         CC(2,1,i)= b(9,i)
 966         CC(1,2,i)= b(9,i)
 967         CC(1,1,-i)= b(7,i)
 968         CC(2,2,-i)=-b(7,i)
 969         CC(2,1,-i)=-b(9,i)
 970         CC(1,2,-i)=-b(9,i)
 971 c        CC(1,1,i)=0.0d0
 972 c        CC(2,2,i)=0.0d0
 973 c        CC(2,1,i)=0.0d0
 974 c        CC(1,2,i)=0.0d0
 975         Ctilde(1,1,i)=b(7,i)
 976         Ctilde(1,2,i)=b(9,i)
 977         Ctilde(2,1,i)=-b(9,i)
 978         Ctilde(2,2,i)=b(7,i)
 979         Ctilde(1,1,-i)=b(7,i)
 980         Ctilde(1,2,-i)=-b(9,i)
 981         Ctilde(2,1,-i)=b(9,i)
 982         Ctilde(2,2,-i)=b(7,i)
 983
 984 c        Ctilde(1,1,i)=0.0d0
 985 c        Ctilde(1,2,i)=0.0d0
 986 c        Ctilde(2,1,i)=0.0d0
 987 c        Ctilde(2,2,i)=0.0d0
 988         DD(1,1,i)= b(6,i)
 989         DD(2,2,i)=-b(6,i)
 990         DD(2,1,i)= b(8,i)
 991         DD(1,2,i)= b(8,i)
 992         DD(1,1,-i)= b(6,i)
 993         DD(2,2,-i)=-b(6,i)
 994         DD(2,1,-i)=-b(8,i)
 995         DD(1,2,-i)=-b(8,i)
 996 c        DD(1,1,i)=0.0d0
 997 c        DD(2,2,i)=0.0d0
 998 c        DD(2,1,i)=0.0d0
 999 c        DD(1,2,i)=0.0d0
1000         Dtilde(1,1,i)=b(6,i)
1001         Dtilde(1,2,i)=b(8,i)
1002         Dtilde(2,1,i)=-b(8,i)
1003         Dtilde(2,2,i)=b(6,i)
1004         Dtilde(1,1,-i)=b(6,i)
1005         Dtilde(1,2,-i)=-b(8,i)
1006         Dtilde(2,1,-i)=b(8,i)
1007         Dtilde(2,2,-i)=b(6,i)
1008
1009 c        Dtilde(1,1,i)=0.0d0
1010 c        Dtilde(1,2,i)=0.0d0
1011 c        Dtilde(2,1,i)=0.0d0
1012 c        Dtilde(2,2,i)=0.0d0
1013         EEold(1,1,i)= b(10,i)+b(11,i)
1014         EEold(2,2,i)=-b(10,i)+b(11,i)
1015         EEold(2,1,i)= b(12,i)-b(13,i)
1016         EEold(1,2,i)= b(12,i)+b(13,i)
1017         EEold(1,1,-i)= b(10,i)+b(11,i)
1018         EEold(2,2,-i)=-b(10,i)+b(11,i)
1019         EEold(2,1,-i)=-b(12,i)+b(13,i)
1020         EEold(1,2,-i)=-b(12,i)-b(13,i)
1021         write(iout,*) "TU DOCHODZE"
1022         print *,"JESTEM"
1023 c        ee(1,1,i)=1.0d0
1024 c        ee(2,2,i)=1.0d0
1025 c        ee(2,1,i)=0.0d0
1026 c        ee(1,2,i)=0.0d0
1027 c        ee(2,1,i)=ee(1,2,i)
1028       enddo
1029 c      lprint=.true.
1030       if (lprint) then
1031       do i=1,nloctyp
1032         write (iout,*) 'Type',i
1033         write (iout,*) 'B1'
1034         write(iout,*) B1(1,i),B1(2,i)
1035         write (iout,*) 'B2'
1036         write(iout,*) B2(1,i),B2(2,i)
1037         write (iout,*) 'CC'
1038         do j=1,2
1039           write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
1040         enddo
1041         write(iout,*) 'DD'
1042         do j=1,2
1043           write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
1044         enddo
1045         write(iout,*) 'EE'
1046         do j=1,2
1047           write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i)
1048         enddo
1049       enddo
1050       endif
1051 c      lprint=.false.
1052
1053 C
1054 C Read electrostatic-interaction parameters
1055 C
1056       if (lprint) then
1057         write (iout,*)
1058         write (iout,'(/a)') 'Electrostatic interaction constants:'
1059         write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
1060      &            'IT','JT','APP','BPP','AEL6','AEL3'
1061       endif
1062       read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
1063       read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
1064       read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
1065       read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
1066       close (ielep)
1067       do i=1,2
1068         do j=1,2
1069         rri=rpp(i,j)**6
1070         app (i,j)=epp(i,j)*rri*rri
1071         bpp (i,j)=-2.0D0*epp(i,j)*rri
1072         ael6(i,j)=elpp6(i,j)*4.2D0**6
1073         ael3(i,j)=elpp3(i,j)*4.2D0**3
1074 c        lprint=.true.
1075         if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
1076      &                    ael6(i,j),ael3(i,j)
1077 c        lprint=.false.
1078         enddo
1079       enddo
1080 C
1081 C Read side-chain interaction parameters.
1082 C
1083       read (isidep,*,end=117,err=117) ipot,expon
1084       if (ipot.lt.1 .or. ipot.gt.5) then
1085         write (iout,'(2a)') 'Error while reading SC interaction',
1086      &               'potential file - unknown potential type.'
1087 #ifdef MPI
1088         call MPI_Finalize(Ierror)
1089 #endif
1090         stop
1091       endif
1092       expon2=expon/2
1093       if(me.eq.king)
1094      & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
1095      & ', exponents are ',expon,2*expon
1096       goto (10,20,30,30,40) ipot
1097 C----------------------- LJ potential ---------------------------------
1098    10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1099      &   (sigma0(i),i=1,ntyp)
1100       if (lprint) then
1101         write (iout,'(/a/)') 'Parameters of the LJ potential:'
1102         write (iout,'(a/)') 'The epsilon array:'
1103         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1104         write (iout,'(/a)') 'One-body parameters:'
1105         write (iout,'(a,4x,a)') 'residue','sigma'
1106         write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1107       endif
1108       goto 50
1109 C----------------------- LJK potential --------------------------------
1110    20 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1111      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1112       if (lprint) then
1113         write (iout,'(/a/)') 'Parameters of the LJK potential:'
1114         write (iout,'(a/)') 'The epsilon array:'
1115         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1116         write (iout,'(/a)') 'One-body parameters:'
1117         write (iout,'(a,4x,2a)') 'residue','   sigma  ','    r0    '
1118         write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
1119      &        i=1,ntyp)
1120       endif
1121       goto 50
1122 C---------------------- GB or BP potential -----------------------------
1123    30 do i=1,ntyp
1124        read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
1125       enddo
1126       read (isidep,*,end=116,err=116)(sigma0(i),i=1,ntyp)
1127       read (isidep,*,end=116,err=116)(sigii(i),i=1,ntyp)
1128       read (isidep,*,end=116,err=116)(chip(i),i=1,ntyp)
1129       read (isidep,*,end=116,err=116)(alp(i),i=1,ntyp)
1130 C now we start reading lipid
1131       do i=1,ntyp
1132        read (isidep,*,end=1161,err=1161)(epslip(i,j),j=i,ntyp)
1133        print *,"WARNING!!"
1134        do j=1,ntyp
1135        epslip(i,j)=epslip(i,j)+0.05d0
1136        enddo
1137       enddo
1138 C For the GB potential convert sigma'**2 into chi'
1139       if (ipot.eq.4) then
1140         do i=1,ntyp
1141           chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1142         enddo
1143       endif
1144       if (lprint) then
1145         write (iout,'(/a/)') 'Parameters of the BP potential:'
1146         write (iout,'(a/)') 'The epsilon array:'
1147         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1148         write (iout,'(/a)') 'One-body parameters:'
1149         write (iout,'(a,4x,4a)') 'residue','   sigma  ','s||/s_|_^2',
1150      &       '    chip  ','    alph  '
1151         write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1152      &                     chip(i),alp(i),i=1,ntyp)
1153       endif
1154       goto 50
1155 C--------------------- GBV potential -----------------------------------
1156    40 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1157      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1158      &  (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1159       if (lprint) then
1160         write (iout,'(/a/)') 'Parameters of the GBV potential:'
1161         write (iout,'(a/)') 'The epsilon array:'
1162         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1163         write (iout,'(/a)') 'One-body parameters:'
1164         write (iout,'(a,4x,5a)') 'residue','   sigma  ','    r0    ',
1165      &      's||/s_|_^2','    chip  ','    alph  '
1166         write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1167      &           sigii(i),chip(i),alp(i),i=1,ntyp)
1168       endif
1169    50 continue
1170       close (isidep)
1171 C-----------------------------------------------------------------------
1172 C Calculate the "working" parameters of SC interactions.
1173       do i=2,ntyp
1174         do j=1,i-1
1175           eps(i,j)=eps(j,i)
1176           epslip(i,j)=epslip(j,i)
1177         enddo
1178       enddo
1179       do i=1,ntyp
1180         do j=i,ntyp
1181           sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1182           sigma(j,i)=sigma(i,j)
1183           rs0(i,j)=dwa16*sigma(i,j)
1184           rs0(j,i)=rs0(i,j)
1185         enddo
1186       enddo
1187       if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1188      & 'Working parameters of the SC interactions:',
1189      & '     a    ','     b    ','   augm   ','  sigma ','   r0   ',
1190      & '  chi1   ','   chi2   '
1191       do i=1,ntyp
1192         do j=i,ntyp
1193           epsij=eps(i,j)
1194           if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1195             rrij=sigma(i,j)
1196           else
1197             rrij=rr0(i)+rr0(j)
1198           endif
1199           r0(i,j)=rrij
1200           r0(j,i)=rrij
1201           rrij=rrij**expon
1202           epsij=eps(i,j)
1203           sigeps=dsign(1.0D0,epsij)
1204           epsij=dabs(epsij)
1205           aa_aq(i,j)=epsij*rrij*rrij
1206           bb_aq(i,j)=-sigeps*epsij*rrij
1207           aa_aq(j,i)=aa_aq(i,j)
1208           bb_aq(j,i)=bb_aq(i,j)
1209           epsijlip=epslip(i,j)
1210           sigeps=dsign(1.0D0,epsijlip)
1211           epsijlip=dabs(epsijlip)
1212           aa_lip(i,j)=epsijlip*rrij*rrij
1213           bb_lip(i,j)=-sigeps*epsijlip*rrij
1214           aa_lip(j,i)=aa_lip(i,j)
1215           bb_lip(j,i)=bb_lip(i,j)
1216           if (ipot.gt.2) then
1217             sigt1sq=sigma0(i)**2
1218             sigt2sq=sigma0(j)**2
1219             sigii1=sigii(i)
1220             sigii2=sigii(j)
1221             ratsig1=sigt2sq/sigt1sq
1222             ratsig2=1.0D0/ratsig1
1223             chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1224             if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1225             rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1226           else
1227             rsum_max=sigma(i,j)
1228           endif
1229 c         if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1230             sigmaii(i,j)=rsum_max
1231             sigmaii(j,i)=rsum_max
1232 c         else
1233 c           sigmaii(i,j)=r0(i,j)
1234 c           sigmaii(j,i)=r0(i,j)
1235 c         endif
1236 cd        write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1237           if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1238             r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1239             augm(i,j)=epsij*r_augm**(2*expon)
1240 c           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1241             augm(j,i)=augm(i,j)
1242           else
1243             augm(i,j)=0.0D0
1244             augm(j,i)=0.0D0
1245           endif
1246           if (lprint) then
1247             write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1248      &      restyp(i),restyp(j),aa_aq(i,j),bb_aq(i,j),augm(i,j),
1249      &      sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1250           endif
1251         enddo
1252       enddo
1253 #ifdef OLDSCP
1254 C
1255 C Define the SC-p interaction constants (hard-coded; old style)
1256 C
1257       do i=1,ntyp
1258 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1259 C helix formation)
1260 c       aad(i,1)=0.3D0*4.0D0**12
1261 C Following line for constants currently implemented
1262 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1263         aad(i,1)=1.5D0*4.0D0**12
1264 c       aad(i,1)=0.17D0*5.6D0**12
1265         aad(i,2)=aad(i,1)
1266 C "Soft" SC-p repulsion
1267         bad(i,1)=0.0D0
1268 C Following line for constants currently implemented
1269 c       aad(i,1)=0.3D0*4.0D0**6
1270 C "Hard" SC-p repulsion
1271         bad(i,1)=3.0D0*4.0D0**6
1272 c       bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1273         bad(i,2)=bad(i,1)
1274 c       aad(i,1)=0.0D0
1275 c       aad(i,2)=0.0D0
1276 c       bad(i,1)=1228.8D0
1277 c       bad(i,2)=1228.8D0
1278       enddo
1279 #else
1280 C
1281 C 8/9/01 Read the SC-p interaction constants from file
1282 C
1283       do i=1,ntyp
1284         read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1285       enddo
1286       do i=1,ntyp
1287         aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1288         aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1289         bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1290         bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1291       enddo
1292 c      lprint=.true.
1293       if (lprint) then
1294         write (iout,*) "Parameters of SC-p interactions:"
1295         do i=1,ntyp
1296           write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1297      &     eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1298         enddo
1299       endif
1300 c      lprint=.false.
1301 #endif
1302 C
1303 C Define the constants of the disulfide bridge
1304 C
1305 C      ebr=-12.00D0
1306 c
1307 c Old arbitrary potential - commented out.
1308 c
1309 c      dbr= 4.20D0
1310 c      fbr= 3.30D0
1311 c
1312 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1313 c energy surface of diethyl disulfide.
1314 c A. Liwo and U. Kozlowska, 11/24/03
1315 c
1316 C      D0CM = 3.78d0
1317 C      AKCM = 15.1d0
1318 C      AKTH = 11.0d0
1319 C      AKCT = 12.0d0
1320 C      V1SS =-1.08d0
1321 C      V2SS = 7.61d0
1322 C      V3SS = 13.7d0
1323 c      akcm=0.0d0
1324 c      akth=0.0d0
1325 c      akct=0.0d0
1326 c      v1ss=0.0d0
1327 c      v2ss=0.0d0
1328 c      v3ss=0.0d0
1329
1330 C      if(me.eq.king) then
1331 C      write (iout,'(/a)') "Disulfide bridge parameters:"
1332 C      write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1333 C      write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1334 C      write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1335 C      write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1336 C     &  ' v3ss:',v3ss
1337 C      endif
1338 C set the variables used for shielding effect
1339 C      write (iout,*) "SHIELD MODE",shield_mode
1340 C      if (shield_mode.gt.0) then
1341 C VSolvSphere the volume of solving sphere
1342 C      print *,pi,"pi"
1343 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
1344 C there will be no distinction between proline peptide group and normal peptide
1345 C group in case of shielding parameters
1346 C      VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
1347 C      VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
1348 C      write (iout,*) VSolvSphere,VSolvSphere_div
1349 C long axis of side chain
1350 C      do i=1,ntyp
1351 C      long_r_sidechain(i)=vbldsc0(1,i)
1352 C      short_r_sidechain(i)=sigma0(i)
1353 C      enddo
1354 C lets set the buffor value
1355 C      buff_shield=1.0d0
1356 C      endif
1357       return
1358   111 write (iout,*) "Error reading bending energy parameters."
1359       goto 999
1360   112 write (iout,*) "Error reading rotamer energy parameters."
1361       goto 999
1362   113 write (iout,*) "Error reading torsional energy parameters."
1363       goto 999
1364   114 write (iout,*) "Error reading double torsional energy parameters."
1365       goto 999
1366   115 write (iout,*)
1367      &  "Error reading cumulant (multibody energy) parameters."
1368       goto 999
1369   116 write (iout,*) "Error reading electrostatic energy parameters."
1370       goto 999
1371  1161 write (iout,*) "Error reading electrostatic energy parameters.Lip"
1372       goto 999
1373   117 write (iout,*) "Error reading side chain interaction parameters."
1374       goto 999
1375   118 write (iout,*) "Error reading SCp interaction parameters."
1376       goto 999
1377   119 write (iout,*) "Error reading SCCOR parameters"
1378   999 continue
1379 #ifdef MPI
1380       call MPI_Finalize(Ierror)
1381 #endif
1382       stop
1383       return
1384       end
1385
1386
1387       subroutine getenv_loc(var, val)
1388       character(*) var, val
1389
1390 #ifdef WINIFL
1391       character(2000) line
1392       external ilen
1393
1394       open (196,file='env',status='old',readonly,shared)
1395       iread=0
1396 c      write(*,*)'looking for ',var
1397 10    read(196,*,err=11,end=11)line
1398       iread=index(line,var)
1399 c      write(*,*)iread,' ',var,' ',line
1400       if (iread.eq.0) go to 10
1401 c      write(*,*)'---> ',line
1402 11    continue
1403       if(iread.eq.0) then
1404 c       write(*,*)'CHUJ'
1405        val=''
1406       else
1407        iread=iread+ilen(var)+1
1408        read (line(iread:),*,err=12,end=12) val
1409 c       write(*,*)'OK: ',var,' = ',val
1410       endif
1411       close(196)
1412       return
1413 12    val=''
1414       close(196)
1415 #elif (defined CRAY)
1416       integer lennam,lenval,ierror
1417 c
1418 c        getenv using a POSIX call, useful on the T3D
1419 c        Sept 1996, comment out error check on advice of H. Pritchard
1420 c
1421       lennam = len(var)
1422       if(lennam.le.0) stop '--error calling getenv--'
1423       call pxfgetenv(var,lennam,val,lenval,ierror)
1424 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1425 #else
1426       call getenv(var,val)
1427 #endif
1428 C set the variables used for shielding effect
1429 C      if (shield_mode.gt.0) then
1430 C VSolvSphere the volume of solving sphere
1431 C      print *,pi,"pi"
1432 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
1433 C there will be no distinction between proline peptide group and normal peptide
1434 C group in case of shielding parameters
1435 C      VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
1436 C      VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
1437 C long axis of side chain
1438 C      do i=1,ntyp
1439 C      long_r_sidechain(i)=vbldsc0(1,i)
1440 C      short_r_sidechain(i)=sigma0(i)
1441 C      enddo
1442 C lets set the buffor value
1443 C      buff_shield=1.0d0
1444 C      endif
1445       return
1446       end