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