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       itortyp(ntyp1)=ntortyp
 634       do iblock=1,2
 635       do i=-ntyp,-1
 636        itortyp(i)=-itortyp(-i)
 637       enddo
 638       write (iout,*) 'ntortyp',ntortyp
 639       do i=0,ntortyp-1
 640         do j=-ntortyp+1,ntortyp-1
 641           read (itorp,*,end=113,err=113) nterm(i,j,iblock),
 642      &          nlor(i,j,iblock)
 643           nterm(-i,-j,iblock)=nterm(i,j,iblock)
 644           nlor(-i,-j,iblock)=nlor(i,j,iblock)
 645           v0ij=0.0d0
 646           si=-1.0d0
 647           do k=1,nterm(i,j,iblock)
 648             read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),
 649      &      v2(k,i,j,iblock)
 650             v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
 651             v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
 652             v0ij=v0ij+si*v1(k,i,j,iblock)
 653             si=-si
 654 c         write(iout,*) i,j,k,iblock,nterm(i,j,iblock)
 655 c         write(iout,*) v1(k,-i,-j,iblock),v1(k,i,j,iblock),
 656 c      &v2(k,-i,-j,iblock),v2(k,i,j,iblock)
 657           enddo
 658           do k=1,nlor(i,j,iblock)
 659             read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),
 660      &        vlor2(k,i,j),vlor3(k,i,j)
 661             v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
 662           enddo
 663           v0(i,j,iblock)=v0ij
 664           v0(-i,-j,iblock)=v0ij
 665         enddo
 666       enddo
 667       enddo
 668       close (itorp)
 669       if (lprint) then
 670         write (iout,'(/a/)') 'Torsional constants:'
 671         do iblock=1,2
 672         do i=0,ntortyp-1
 673           do j=-ntortyp+1,ntortyp-1
 674             write (iout,*) 'ityp',i,' jtyp',j
 675             write (iout,*) 'Fourier constants'
 676             do k=1,nterm(i,j,iblock)
 677               write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
 678      &        v2(k,i,j,iblock)
 679             enddo
 680             write (iout,*) 'Lorenz constants'
 681             do k=1,nlor(i,j,iblock)
 682               write (iout,'(3(1pe15.5))')
 683      &         vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
 684             enddo
 685           enddo
 686         enddo
 687         enddo
 688       endif
 689
 690 C
 691 C 6/23/01 Read parameters for double torsionals
 692 C
 693       do iblock=1,2
 694       do i=0,ntortyp-1
 695         do j=-ntortyp+1,ntortyp-1
 696           do k=-ntortyp+1,ntortyp-1
 697             read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
 698 c              write (iout,*) "OK onelett",
 699 c     &         i,j,k,t1,t2,t3
 700
 701             if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
 702      &        .or. t3.ne.toronelet(k)) then
 703               write (iout,*) "Error in double torsional parameter file",
 704      &         i,j,k,t1,t2,t3
 705 #ifdef MPI
 706               call MPI_Finalize(Ierror)
 707 #endif
 708                stop "Error in double torsional parameter file"
 709             endif
 710            read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),
 711      &         ntermd_2(i,j,k,iblock)
 712             ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
 713             ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
 714             read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,
 715      &         ntermd_1(i,j,k,iblock))
 716             read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,
 717      &         ntermd_1(i,j,k,iblock))
 718             read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,
 719      &         ntermd_1(i,j,k,iblock))
 720             read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,
 721      &         ntermd_1(i,j,k,iblock))
 722 C Martix of D parameters for one dimesional foureir series
 723             do l=1,ntermd_1(i,j,k,iblock)
 724              v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
 725              v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
 726              v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
 727              v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
 728 c            write(iout,*) "whcodze" ,
 729 c     & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
 730             enddo
 731             read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),
 732      &         v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
 733      &         v2s(m,l,i,j,k,iblock),
 734      &         m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
 735 C Martix of D parameters for two dimesional fourier series
 736             do l=1,ntermd_2(i,j,k,iblock)
 737              do m=1,l-1
 738              v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
 739              v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
 740              v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
 741              v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
 742              enddo!m
 743             enddo!l
 744           enddo!k
 745         enddo!j
 746       enddo!i
 747       enddo!iblock
 748       if (lprint) then
 749       write (iout,*)
 750       write (iout,*) 'Constants for double torsionals'
 751       do iblock=1,2
 752       do i=0,ntortyp-1
 753         do j=-ntortyp+1,ntortyp-1
 754           do k=-ntortyp+1,ntortyp-1
 755             write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
 756      &        ' nsingle',ntermd_1(i,j,k,iblock),
 757      &        ' ndouble',ntermd_2(i,j,k,iblock)
 758             write (iout,*)
 759             write (iout,*) 'Single angles:'
 760             do l=1,ntermd_1(i,j,k,iblock)
 761               write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
 762      &           v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
 763      &           v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),
 764      &           v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
 765             enddo
 766             write (iout,*)
 767             write (iout,*) 'Pairs of angles:'
 768             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
 769             do l=1,ntermd_2(i,j,k,iblock)
 770               write (iout,'(i5,20f10.5)')
 771      &         l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
 772             enddo
 773             write (iout,*)
 774             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
 775             do l=1,ntermd_2(i,j,k,iblock)
 776               write (iout,'(i5,20f10.5)')
 777      &         l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),
 778      &         (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
 779             enddo
 780             write (iout,*)
 781           enddo
 782         enddo
 783       enddo
 784       enddo
 785       endif
 786 #endif
 787 C read Czybyshev torsional parameters
 788       read (itorkcc,*,end=121,err=121) nkcctyp
 789       read (itorkcc,*,end=121,err=121) (itortyp_kcc(i),i=1,ntyp)
 790       do i=-ntyp,-1
 791         itortyp_kcc(i)=-itortyp_kcc(-i)
 792       enddo
 793       do i=0,nkcctyp
 794         do j=0,nkcctyp
 795 C first we read the cos and sin gamma parameters
 796           read (itorkcc,*,end=121,err=121)
 797      &    nterm_kcc(j,i),nterm_kcc_Tb(j,i)
 798 C           read (itorkcc,*,end=121,err=121) nterm_kcc_Tb(j,i)
 799           do k=1,nterm_kcc(j,i)
 800             do l=1,nterm_kcc_Tb(j,i)
 801               read (itorkcc,*,end=121,err=121) v11_chyb(l,k,j,i)
 802             enddo
 803             do l=1,nterm_kcc_Tb(j,i)
 804               read (itorkcc,*,end=121,err=121) v21_chyb(l,k,j,i)
 805             enddo
 806             do l=1,nterm_kcc_Tb(j,i)
 807               read (itorkcc,*,end=121,err=121) v12_chyb(l,k,j,i)
 808             enddo
 809             do l=1,nterm_kcc_Tb(j,i)
 810               read (itorkcc,*,end=121,err=121) v22_chyb(l,k,j,i)
 811             enddo
 812             read (itorkcc,*,end=121,err=121) v1_kcc(k,j,i)
 813             read (itorkcc,*,end=121,err=121) v2_kcc(k,j,i)
 814           enddo
 815         enddo
 816       enddo
 817       if (lprint) then
 818 c Print valence-torsional parameters
 819         write (iout,'(a)')
 820      &    "Parameters of the valence-torsional potentials"
 821         do i=0,nkcctyp
 822         do j=0,nkcctyp
 823         write (iout,'(3a)') "Type ",toronelet(i),toronelet(j)
 824         write (iout,'(2a20,a15)') "v_kcc","v1_chyb","v2_chyb"
 825         do k=1,nterm_kcc(j,i)
 826           write (iout,'(i5,f15.10,i5,2f15.10)')
 827      &      k,v1_kcc(k,j,i),1,v11_chyb(1,k,j,i),v21_chyb(1,k,j,i)
 828           do l=2,nterm_kcc_Tb(j,i)
 829             write (iout,'(20x,i5,2f15.10)')
 830      &        l,v11_chyb(l,k,j,i),v21_chyb(l,k,j,i)
 831           enddo
 832           write (iout,'(i5,f15.10,i5,2f15.10)')
 833      &      k,v2_kcc(k,j,i),1,v12_chyb(1,k,j,i),v22_chyb(1,k,j,i)
 834           do l=2,nterm_kcc_Tb(j,i)
 835             write (iout,'(20x,i5,2f15.10)')
 836      &        l,v12_chyb(l,k,j,i),v22_chyb(l,k,j,i)
 837           enddo
 838           write (iout,'(a)')
 839         enddo
 840         enddo
 841         enddo
 842       endif
 843 C here will be the apropriate recalibrating for D-aminoacid
 844 C        read (ithetkcc,*,end=121,err=121) nkcctyp
 845       do i=0,nkcctyp
 846         read (ithetkcc,*,end=121,err=121) nbend_kcc_Tb(i)
 847         do j=1,nbend_kcc_Tb(i)
 848           read (ithetkcc,*,end=121,err=121) v1bend_chyb(j,i)
 849         enddo
 850       enddo
 851       if (lprint) then
 852         write (iout,'(a)')
 853      &    "Parameters of the valence-only potentials"
 854         do i=0,nkcctyp
 855         write (iout,'(2a)') "Type ",toronelet(i)
 856         do k=1,nbend_kcc_Tb(i)
 857           write(iout,'(i5,f15.10)') k,v1bend_chyb(k,i)
 858         enddo
 859         enddo
 860       endif
 861 C Read of Side-chain backbone correlation parameters
 862 C Modified 11 May 2012 by Adasko
 863 CCC
 864 C
 865       read (isccor,*,end=119,err=119) nsccortyp
 866 #ifdef SCCORPDB
 867       read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
 868       do i=-ntyp,-1
 869         isccortyp(i)=-isccortyp(-i)
 870       enddo
 871       iscprol=isccortyp(20)
 872 c      write (iout,*) 'ntortyp',ntortyp
 873       maxinter=3
 874 cc maxinter is maximum interaction sites
 875       do l=1,maxinter
 876       do i=1,nsccortyp
 877         do j=1,nsccortyp
 878           read (isccor,*,end=119,err=119)
 879      &nterm_sccor(i,j),nlor_sccor(i,j)
 880           v0ijsccor=0.0d0
 881           v0ijsccor1=0.0d0
 882           v0ijsccor2=0.0d0
 883           v0ijsccor3=0.0d0
 884           si=-1.0d0
 885           nterm_sccor(-i,j)=nterm_sccor(i,j)
 886           nterm_sccor(-i,-j)=nterm_sccor(i,j)
 887           nterm_sccor(i,-j)=nterm_sccor(i,j)
 888           do k=1,nterm_sccor(i,j)
 889             read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
 890      &    ,v2sccor(k,l,i,j)
 891             if (j.eq.iscprol) then
 892              if (i.eq.isccortyp(10)) then
 893              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
 894              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
 895              else
 896              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
 897      &                        +v2sccor(k,l,i,j)*dsqrt(0.75d0)
 898              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
 899      &                        +v1sccor(k,l,i,j)*dsqrt(0.75d0)
 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             else
 906              if (i.eq.isccortyp(10)) then
 907              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
 908              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
 909              else
 910                if (j.eq.isccortyp(10)) then
 911              v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
 912              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
 913                else
 914              v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
 915              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
 916              v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
 917              v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
 918              v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
 919              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
 920                 endif
 921                endif
 922             endif
 923             v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
 924             v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
 925             v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
 926             v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
 927             si=-si
 928           enddo
 929           do k=1,nlor_sccor(i,j)
 930             read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
 931      &        vlor2sccor(k,i,j),vlor3sccor(k,i,j)
 932             v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
 933      &(1+vlor3sccor(k,i,j)**2)
 934           enddo
 935           v0sccor(l,i,j)=v0ijsccor
 936           v0sccor(l,-i,j)=v0ijsccor1
 937           v0sccor(l,i,-j)=v0ijsccor2
 938           v0sccor(l,-i,-j)=v0ijsccor3
 939         enddo
 940       enddo
 941       enddo
 942       close (isccor)
 943 #else
 944       read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
 945 c      write (iout,*) 'ntortyp',ntortyp
 946       maxinter=3
 947 cc maxinter is maximum interaction sites
 948       do l=1,maxinter
 949       do i=1,nsccortyp
 950         do j=1,nsccortyp
 951           read (isccor,*,end=119,err=119)
 952      & nterm_sccor(i,j),nlor_sccor(i,j)
 953           v0ijsccor=0.0d0
 954           si=-1.0d0
 955
 956           do k=1,nterm_sccor(i,j)
 957             read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
 958      &    ,v2sccor(k,l,i,j)
 959             v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
 960             si=-si
 961           enddo
 962           do k=1,nlor_sccor(i,j)
 963             read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
 964      &        vlor2sccor(k,i,j),vlor3sccor(k,i,j)
 965             v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
 966      &(1+vlor3sccor(k,i,j)**2)
 967           enddo
 968           v0sccor(l,i,j)=v0ijsccor
 969         enddo
 970       enddo
 971       enddo
 972       close (isccor)
 973
 974 #endif
 975       if (lprint) then
 976         write (iout,'(/a/)') 'Torsional constants:'
 977         do l=1,maxinter
 978         do i=1,nsccortyp
 979           do j=1,nsccortyp
 980             write (iout,*) 'ityp',i,' jtyp',j
 981             write (iout,*) 'Fourier constants'
 982             do k=1,nterm_sccor(i,j)
 983       write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
 984             enddo
 985             write (iout,*) 'Lorenz constants'
 986             do k=1,nlor_sccor(i,j)
 987               write (iout,'(3(1pe15.5))')
 988      &         vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
 989             enddo
 990           enddo
 991         enddo
 992         enddo
 993       endif
 994
 995 C
 996 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
 997 C         interaction energy of the Gly, Ala, and Pro prototypes.
 998 C
 999       if (lprint) then
1000         write (iout,*)
1001         write (iout,*) "Coefficients of the cumulants"
1002       endif
1003       read (ifourier,*) nloctyp
1004 #ifdef NEWCORR
1005       read (ifourier,*,end=115,err=115) (itype2loc(i),i=1,ntyp)
1006       read (ifourier,*,end=115,err=115) (iloctyp(i),i=0,nloctyp-1)
1007       itype2loc(ntyp1)=nloctyp
1008       iloctyp(nloctyp)=ntyp1
1009 #else
1010       do i=1,ntyp1
1011         itype2loc(i)=itortyp(i)
1012       enddo
1013       iloctyp(0)=10
1014       iloctyp(1)=9
1015       iloctyp(2)=20
1016       iloctyp(3)=ntyp1
1017 #endif
1018       do i=1,ntyp1
1019         itype2loc(-i)=-itype2loc(i)
1020       enddo
1021       do i=1,nloctyp
1022         iloctyp(-i)=-iloctyp(i)
1023       enddo
1024       write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1)
1025       write (iout,*) "nloctyp",nloctyp,
1026      &  " iloctyp",(iloctyp(i),i=0,nloctyp)
1027       do i=0,nloctyp-1
1028         read (ifourier,*,end=115,err=115)
1029         read (ifourier,*,end=115,err=115) (b(ii,i),ii=1,13)
1030 #ifdef NEWCORR
1031         read (ifourier,*,end=115,err=115) (bnew1(ii,1,i),ii=1,3)
1032         read (ifourier,*,end=115,err=115) (bnew2(ii,1,i),ii=1,3)
1033         read (ifourier,*,end=115,err=115) (bnew1(ii,2,i),ii=1,1)
1034         read (ifourier,*,end=115,err=115) (bnew2(ii,2,i),ii=1,1)
1035         read (ifourier,*,end=115,err=115) (eenew(ii,i),ii=1,1)
1036 #endif
1037         if (lprint) then
1038         write (iout,*) 'Type',i
1039         write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
1040         endif
1041 c        B1(1,i)  = b(3)
1042 c        B1(2,i)  = b(5)
1043 c        B1(1,-i) = b(3)
1044 c        B1(2,-i) = -b(5)
1045 c        b1(1,i)=0.0d0
1046 c        b1(2,i)=0.0d0
1047 c        B1tilde(1,i) = b(3)
1048 c        B1tilde(2,i) =-b(5)
1049 c        B1tilde(1,-i) =-b(3)
1050 c        B1tilde(2,-i) =b(5)
1051 c        b1tilde(1,i)=0.0d0
1052 c        b1tilde(2,i)=0.0d0
1053 c        B2(1,i)  = b(2)
1054 c        B2(2,i)  = b(4)
1055 c        B2(1,-i)  =b(2)
1056 c        B2(2,-i)  =-b(4)
1057 cc        B1tilde(1,i) = b(3,i)
1058 cc        B1tilde(2,i) =-b(5,i)
1059 C        B1tilde(1,-i) =-b(3,i)
1060 C        B1tilde(2,-i) =b(5,i)
1061 cc        b1tilde(1,i)=0.0d0
1062 cc        b1tilde(2,i)=0.0d0
1063 cc        B2(1,i)  = b(2,i)
1064 cc        B2(2,i)  = b(4,i)
1065 C        B2(1,-i)  =b(2,i)
1066 C        B2(2,-i)  =-b(4,i)
1067
1068 c        b2(1,i)=0.0d0
1069 c        b2(2,i)=0.0d0
1070         CC(1,1,i)= b(7,i)
1071         CC(2,2,i)=-b(7,i)
1072         CC(2,1,i)= b(9,i)
1073         CC(1,2,i)= b(9,i)
1074         CC(1,1,-i)= b(7,i)
1075         CC(2,2,-i)=-b(7,i)
1076         CC(2,1,-i)=-b(9,i)
1077         CC(1,2,-i)=-b(9,i)
1078 c        CC(1,1,i)=0.0d0
1079 c        CC(2,2,i)=0.0d0
1080 c        CC(2,1,i)=0.0d0
1081 c        CC(1,2,i)=0.0d0
1082         Ctilde(1,1,i)=b(7,i)
1083         Ctilde(1,2,i)=b(9,i)
1084         Ctilde(2,1,i)=-b(9,i)
1085         Ctilde(2,2,i)=b(7,i)
1086         Ctilde(1,1,-i)=b(7,i)
1087         Ctilde(1,2,-i)=-b(9,i)
1088         Ctilde(2,1,-i)=b(9,i)
1089         Ctilde(2,2,-i)=b(7,i)
1090
1091 c        Ctilde(1,1,i)=0.0d0
1092 c        Ctilde(1,2,i)=0.0d0
1093 c        Ctilde(2,1,i)=0.0d0
1094 c        Ctilde(2,2,i)=0.0d0
1095         DD(1,1,i)= b(6,i)
1096         DD(2,2,i)=-b(6,i)
1097         DD(2,1,i)= b(8,i)
1098         DD(1,2,i)= b(8,i)
1099         DD(1,1,-i)= b(6,i)
1100         DD(2,2,-i)=-b(6,i)
1101         DD(2,1,-i)=-b(8,i)
1102         DD(1,2,-i)=-b(8,i)
1103 c        DD(1,1,i)=0.0d0
1104 c        DD(2,2,i)=0.0d0
1105 c        DD(2,1,i)=0.0d0
1106 c        DD(1,2,i)=0.0d0
1107         Dtilde(1,1,i)=b(6,i)
1108         Dtilde(1,2,i)=b(8,i)
1109         Dtilde(2,1,i)=-b(8,i)
1110         Dtilde(2,2,i)=b(6,i)
1111         Dtilde(1,1,-i)=b(6,i)
1112         Dtilde(1,2,-i)=-b(8,i)
1113         Dtilde(2,1,-i)=b(8,i)
1114         Dtilde(2,2,-i)=b(6,i)
1115
1116 c        Dtilde(1,1,i)=0.0d0
1117 c        Dtilde(1,2,i)=0.0d0
1118 c        Dtilde(2,1,i)=0.0d0
1119 c        Dtilde(2,2,i)=0.0d0
1120         EEold(1,1,i)= b(10,i)+b(11,i)
1121         EEold(2,2,i)=-b(10,i)+b(11,i)
1122         EEold(2,1,i)= b(12,i)-b(13,i)
1123         EEold(1,2,i)= b(12,i)+b(13,i)
1124         EEold(1,1,-i)= b(10,i)+b(11,i)
1125         EEold(2,2,-i)=-b(10,i)+b(11,i)
1126         EEold(2,1,-i)=-b(12,i)+b(13,i)
1127         EEold(1,2,-i)=-b(12,i)-b(13,i)
1128         write(iout,*) "TU DOCHODZE"
1129         print *,"JESTEM"
1130 c        ee(1,1,i)=1.0d0
1131 c        ee(2,2,i)=1.0d0
1132 c        ee(2,1,i)=0.0d0
1133 c        ee(1,2,i)=0.0d0
1134 c        ee(2,1,i)=ee(1,2,i)
1135       enddo
1136       lprint=.true.
1137       if (lprint) then
1138       do i=1,nloctyp
1139         write (iout,*) 'Type',i
1140         write (iout,*) 'B1'
1141         write(iout,*) B1(1,i),B1(2,i)
1142         write (iout,*) 'B2'
1143         write(iout,*) B2(1,i),B2(2,i)
1144         write (iout,*) 'CC'
1145         do j=1,2
1146           write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
1147         enddo
1148         write(iout,*) 'DD'
1149         do j=1,2
1150           write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
1151         enddo
1152         write(iout,*) 'EE'
1153         do j=1,2
1154           write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i)
1155         enddo
1156       enddo
1157       endif
1158       lprint=.false.
1159
1160 C
1161 C Read electrostatic-interaction parameters
1162 C
1163       if (lprint) then
1164         write (iout,*)
1165         write (iout,'(/a)') 'Electrostatic interaction constants:'
1166         write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
1167      &            'IT','JT','APP','BPP','AEL6','AEL3'
1168       endif
1169       read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
1170       read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
1171       read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
1172       read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
1173       close (ielep)
1174       do i=1,2
1175         do j=1,2
1176         rri=rpp(i,j)**6
1177         app (i,j)=epp(i,j)*rri*rri
1178         bpp (i,j)=-2.0D0*epp(i,j)*rri
1179         ael6(i,j)=elpp6(i,j)*4.2D0**6
1180         ael3(i,j)=elpp3(i,j)*4.2D0**3
1181 c        lprint=.true.
1182         if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
1183      &                    ael6(i,j),ael3(i,j)
1184 c        lprint=.false.
1185         enddo
1186       enddo
1187 C
1188 C Read side-chain interaction parameters.
1189 C
1190       read (isidep,*,end=117,err=117) ipot,expon
1191       if (ipot.lt.1 .or. ipot.gt.5) then
1192         write (iout,'(2a)') 'Error while reading SC interaction',
1193      &               'potential file - unknown potential type.'
1194 #ifdef MPI
1195         call MPI_Finalize(Ierror)
1196 #endif
1197         stop
1198       endif
1199       expon2=expon/2
1200       if(me.eq.king)
1201      & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
1202      & ', exponents are ',expon,2*expon
1203       goto (10,20,30,30,40) ipot
1204 C----------------------- LJ potential ---------------------------------
1205    10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1206      &   (sigma0(i),i=1,ntyp)
1207       if (lprint) then
1208         write (iout,'(/a/)') 'Parameters of the LJ potential:'
1209         write (iout,'(a/)') 'The epsilon array:'
1210         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1211         write (iout,'(/a)') 'One-body parameters:'
1212         write (iout,'(a,4x,a)') 'residue','sigma'
1213         write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1214       endif
1215       goto 50
1216 C----------------------- LJK potential --------------------------------
1217    20 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1218      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1219       if (lprint) then
1220         write (iout,'(/a/)') 'Parameters of the LJK potential:'
1221         write (iout,'(a/)') 'The epsilon array:'
1222         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1223         write (iout,'(/a)') 'One-body parameters:'
1224         write (iout,'(a,4x,2a)') 'residue','   sigma  ','    r0    '
1225         write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
1226      &        i=1,ntyp)
1227       endif
1228       goto 50
1229 C---------------------- GB or BP potential -----------------------------
1230    30 do i=1,ntyp
1231        read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
1232       enddo
1233       read (isidep,*,end=116,err=116)(sigma0(i),i=1,ntyp)
1234       read (isidep,*,end=116,err=116)(sigii(i),i=1,ntyp)
1235       read (isidep,*,end=116,err=116)(chip(i),i=1,ntyp)
1236       read (isidep,*,end=116,err=116)(alp(i),i=1,ntyp)
1237 C now we start reading lipid
1238       do i=1,ntyp
1239        read (isidep,*,end=1161,err=1161)(epslip(i,j),j=i,ntyp)
1240        write(iout,*) "epslip", i, (epslip(i,j),j=i,ntyp)
1241
1242 C       print *,"WARNING!!"
1243 C       do j=1,ntyp
1244 C       epslip(i,j)=epslip(i,j)+0.05d0
1245 C       enddo
1246       enddo
1247       write(iout,*) epslip(1,1),"OK?"
1248 C For the GB potential convert sigma'**2 into chi'
1249       if (ipot.eq.4) then
1250         do i=1,ntyp
1251           chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1252         enddo
1253       endif
1254       if (lprint) then
1255         write (iout,'(/a/)') 'Parameters of the BP potential:'
1256         write (iout,'(a/)') 'The epsilon array:'
1257         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1258         write (iout,'(/a)') 'One-body parameters:'
1259         write (iout,'(a,4x,4a)') 'residue','   sigma  ','s||/s_|_^2',
1260      &       '    chip  ','    alph  '
1261         write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1262      &                     chip(i),alp(i),i=1,ntyp)
1263       endif
1264       goto 50
1265 C--------------------- GBV potential -----------------------------------
1266    40 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1267      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1268      &  (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1269       if (lprint) then
1270         write (iout,'(/a/)') 'Parameters of the GBV potential:'
1271         write (iout,'(a/)') 'The epsilon array:'
1272         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1273         write (iout,'(/a)') 'One-body parameters:'
1274         write (iout,'(a,4x,5a)') 'residue','   sigma  ','    r0    ',
1275      &      's||/s_|_^2','    chip  ','    alph  '
1276         write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1277      &           sigii(i),chip(i),alp(i),i=1,ntyp)
1278       endif
1279    50 continue
1280       close (isidep)
1281 C-----------------------------------------------------------------------
1282 C Calculate the "working" parameters of SC interactions.
1283       do i=2,ntyp
1284         do j=1,i-1
1285           eps(i,j)=eps(j,i)
1286           epslip(i,j)=epslip(j,i)
1287         enddo
1288       enddo
1289       do i=1,ntyp
1290         do j=i,ntyp
1291           sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1292           sigma(j,i)=sigma(i,j)
1293           rs0(i,j)=dwa16*sigma(i,j)
1294           rs0(j,i)=rs0(i,j)
1295         enddo
1296       enddo
1297       if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1298      & 'Working parameters of the SC interactions:',
1299      & '     a    ','     b    ','   augm   ','  sigma ','   r0   ',
1300      & '  chi1   ','   chi2   '
1301       do i=1,ntyp
1302         do j=i,ntyp
1303           epsij=eps(i,j)
1304           if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1305             rrij=sigma(i,j)
1306           else
1307             rrij=rr0(i)+rr0(j)
1308           endif
1309           r0(i,j)=rrij
1310           r0(j,i)=rrij
1311           rrij=rrij**expon
1312           epsij=eps(i,j)
1313           sigeps=dsign(1.0D0,epsij)
1314           epsij=dabs(epsij)
1315           aa_aq(i,j)=epsij*rrij*rrij
1316           bb_aq(i,j)=-sigeps*epsij*rrij
1317           aa_aq(j,i)=aa_aq(i,j)
1318           bb_aq(j,i)=bb_aq(i,j)
1319           epsijlip=epslip(i,j)
1320           sigeps=dsign(1.0D0,epsijlip)
1321           epsijlip=dabs(epsijlip)
1322           aa_lip(i,j)=epsijlip*rrij*rrij
1323           bb_lip(i,j)=-sigeps*epsijlip*rrij
1324           aa_lip(j,i)=aa_lip(i,j)
1325           bb_lip(j,i)=bb_lip(i,j)
1326           if (ipot.gt.2) then
1327             sigt1sq=sigma0(i)**2
1328             sigt2sq=sigma0(j)**2
1329             sigii1=sigii(i)
1330             sigii2=sigii(j)
1331             ratsig1=sigt2sq/sigt1sq
1332             ratsig2=1.0D0/ratsig1
1333             chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1334             if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1335             rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1336           else
1337             rsum_max=sigma(i,j)
1338           endif
1339 c         if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1340             sigmaii(i,j)=rsum_max
1341             sigmaii(j,i)=rsum_max
1342 c         else
1343 c           sigmaii(i,j)=r0(i,j)
1344 c           sigmaii(j,i)=r0(i,j)
1345 c         endif
1346 cd        write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1347           if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1348             r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1349             augm(i,j)=epsij*r_augm**(2*expon)
1350 c           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1351             augm(j,i)=augm(i,j)
1352           else
1353             augm(i,j)=0.0D0
1354             augm(j,i)=0.0D0
1355           endif
1356           if (lprint) then
1357             write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1358      &      restyp(i),restyp(j),aa_aq(i,j),bb_aq(i,j),augm(i,j),
1359      &      sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1360           endif
1361         enddo
1362       enddo
1363       write(iout,*) "tube param"
1364       read(itube,*) epspeptube,sigmapeptube
1365       sigmapeptube=sigmapeptube**6
1366       sigeps=dsign(1.0D0,epspeptube)
1367       epspeptube=dabs(epspeptube)
1368       pep_aa_tube=4.0d0*epspeptube*sigmapeptube**2
1369       pep_bb_tube=-sigeps*4.0d0*epspeptube*sigmapeptube
1370       write(iout,*) pep_aa_tube,pep_bb_tube
1371       do i=1,ntyp
1372        read(itube,*) epssctube,sigmasctube
1373        sigmasctube=sigmasctube**6
1374        sigeps=dsign(1.0D0,epssctube)
1375        epssctube=dabs(epssctube)
1376        sc_aa_tube_par(i)=4.0d0*epssctube*sigmasctube**2
1377        sc_bb_tube_par(i)=-sigeps*4.0d0*epssctube*sigmasctube
1378       write(iout,*) sc_aa_tube_par(i), sc_bb_tube_par(i)
1379       enddo
1380
1381 #ifdef OLDSCP
1382 C
1383 C Define the SC-p interaction constants (hard-coded; old style)
1384 C
1385       do i=1,ntyp
1386 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1387 C helix formation)
1388 c       aad(i,1)=0.3D0*4.0D0**12
1389 C Following line for constants currently implemented
1390 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1391         aad(i,1)=1.5D0*4.0D0**12
1392 c       aad(i,1)=0.17D0*5.6D0**12
1393         aad(i,2)=aad(i,1)
1394 C "Soft" SC-p repulsion
1395         bad(i,1)=0.0D0
1396 C Following line for constants currently implemented
1397 c       aad(i,1)=0.3D0*4.0D0**6
1398 C "Hard" SC-p repulsion
1399         bad(i,1)=3.0D0*4.0D0**6
1400 c       bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1401         bad(i,2)=bad(i,1)
1402 c       aad(i,1)=0.0D0
1403 c       aad(i,2)=0.0D0
1404 c       bad(i,1)=1228.8D0
1405 c       bad(i,2)=1228.8D0
1406       enddo
1407 #else
1408 C
1409 C 8/9/01 Read the SC-p interaction constants from file
1410 C
1411       do i=1,ntyp
1412         read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1413       enddo
1414       do i=1,ntyp
1415         aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1416         aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1417         bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1418         bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1419       enddo
1420 c      lprint=.true.
1421       if (lprint) then
1422         write (iout,*) "Parameters of SC-p interactions:"
1423         do i=1,ntyp
1424           write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1425      &     eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1426         enddo
1427       endif
1428 c      lprint=.false.
1429 #endif
1430 C
1431 C Define the constants of the disulfide bridge
1432 C
1433 C      ebr=-12.00D0
1434 c
1435 c Old arbitrary potential - commented out.
1436 c
1437 c      dbr= 4.20D0
1438 c      fbr= 3.30D0
1439 c
1440 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1441 c energy surface of diethyl disulfide.
1442 c A. Liwo and U. Kozlowska, 11/24/03
1443 c
1444 C      D0CM = 3.78d0
1445 C      AKCM = 15.1d0
1446 C      AKTH = 11.0d0
1447 C      AKCT = 12.0d0
1448 C      V1SS =-1.08d0
1449 C      V2SS = 7.61d0
1450 C      V3SS = 13.7d0
1451 c      akcm=0.0d0
1452 c      akth=0.0d0
1453 c      akct=0.0d0
1454 c      v1ss=0.0d0
1455 c      v2ss=0.0d0
1456 c      v3ss=0.0d0
1457
1458 C      if(me.eq.king) then
1459 C      write (iout,'(/a)') "Disulfide bridge parameters:"
1460 C      write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1461 C      write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1462 C      write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1463 C      write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1464 C     &  ' v3ss:',v3ss
1465 C      endif
1466 C set the variables used for shielding effect
1467 C      write (iout,*) "SHIELD MODE",shield_mode
1468 C      if (shield_mode.gt.0) then
1469 C VSolvSphere the volume of solving sphere
1470 C      print *,pi,"pi"
1471 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
1472 C there will be no distinction between proline peptide group and normal peptide
1473 C group in case of shielding parameters
1474 C      VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
1475 C      VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
1476 C      write (iout,*) VSolvSphere,VSolvSphere_div
1477 C long axis of side chain
1478 C      do i=1,ntyp
1479 C      long_r_sidechain(i)=vbldsc0(1,i)
1480 C      short_r_sidechain(i)=sigma0(i)
1481 C      enddo
1482 C lets set the buffor value
1483 C      buff_shield=1.0d0
1484 C      endif
1485       return
1486   111 write (iout,*) "Error reading bending energy parameters."
1487       goto 999
1488   112 write (iout,*) "Error reading rotamer energy parameters."
1489       goto 999
1490   113 write (iout,*) "Error reading torsional energy parameters."
1491       goto 999
1492   114 write (iout,*) "Error reading double torsional energy parameters."
1493       goto 999
1494   115 write (iout,*)
1495      &  "Error reading cumulant (multibody energy) parameters."
1496       goto 999
1497   116 write (iout,*) "Error reading electrostatic energy parameters."
1498       goto 999
1499  1161 write (iout,*) "Error reading electrostatic energy parameters.Lip"
1500       goto 999
1501   117 write (iout,*) "Error reading side chain interaction parameters."
1502       goto 999
1503   118 write (iout,*) "Error reading SCp interaction parameters."
1504       goto 999
1505   119 write (iout,*) "Error reading SCCOR parameters"
1506       go to 999
1507   121 write (iout,*) "Error in Czybyshev parameters"
1508   999 continue
1509 #ifdef MPI
1510       call MPI_Finalize(Ierror)
1511 #endif
1512       stop
1513       return
1514       end
1515
1516
1517       subroutine getenv_loc(var, val)
1518       character(*) var, val
1519
1520 #ifdef WINIFL
1521       character(2000) line
1522       external ilen
1523
1524       open (196,file='env',status='old',readonly,shared)
1525       iread=0
1526 c      write(*,*)'looking for ',var
1527 10    read(196,*,err=11,end=11)line
1528       iread=index(line,var)
1529 c      write(*,*)iread,' ',var,' ',line
1530       if (iread.eq.0) go to 10
1531 c      write(*,*)'---> ',line
1532 11    continue
1533       if(iread.eq.0) then
1534 c       write(*,*)'CHUJ'
1535        val=''
1536       else
1537        iread=iread+ilen(var)+1
1538        read (line(iread:),*,err=12,end=12) val
1539 c       write(*,*)'OK: ',var,' = ',val
1540       endif
1541       close(196)
1542       return
1543 12    val=''
1544       close(196)
1545 #elif (defined CRAY)
1546       integer lennam,lenval,ierror
1547 c
1548 c        getenv using a POSIX call, useful on the T3D
1549 c        Sept 1996, comment out error check on advice of H. Pritchard
1550 c
1551       lennam = len(var)
1552       if(lennam.le.0) stop '--error calling getenv--'
1553       call pxfgetenv(var,lennam,val,lenval,ierror)
1554 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1555 #else
1556       call getenv(var,val)
1557 #endif
1558 C set the variables used for shielding effect
1559 C      if (shield_mode.gt.0) then
1560 C VSolvSphere the volume of solving sphere
1561 C      print *,pi,"pi"
1562 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
1563 C there will be no distinction between proline peptide group and normal peptide
1564 C group in case of shielding parameters
1565 C      VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
1566 C      VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
1567 C long axis of side chain
1568 C      do i=1,ntyp
1569 C      long_r_sidechain(i)=vbldsc0(1,i)
1570 C      short_r_sidechain(i)=sigma0(i)
1571 C      enddo
1572 C lets set the buffor value
1573 C      buff_shield=1.0d0
1574 C      endif
1575       return
1576       end