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