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