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