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       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       do i=0,nloctyp-1
 917         read (ifourier,*,end=115,err=115)
 918         read (ifourier,*,end=115,err=115) (b(ii),ii=1,13)
 919         if (lprint) then
 920         write (iout,*) 'Type',i
 921         write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
 922         endif
 923         B1(1,i)  = b(3)
 924         B1(2,i)  = b(5)
 925         B1(1,-i) = b(3)
 926         B1(2,-i) = -b(5)
 927 c        b1(1,i)=0.0d0
 928 c        b1(2,i)=0.0d0
 929         B1tilde(1,i) = b(3)
 930         B1tilde(2,i) =-b(5)
 931         B1tilde(1,-i) =-b(3)
 932         B1tilde(2,-i) =b(5)
 933 c        b1tilde(1,i)=0.0d0
 934 c        b1tilde(2,i)=0.0d0
 935         B2(1,i)  = b(2)
 936         B2(2,i)  = b(4)
 937         B2(1,-i)  =b(2)
 938         B2(2,-i)  =-b(4)
 939
 940 c        b2(1,i)=0.0d0
 941 c        b2(2,i)=0.0d0
 942         CC(1,1,i)= b(7)
 943         CC(2,2,i)=-b(7)
 944         CC(2,1,i)= b(9)
 945         CC(1,2,i)= b(9)
 946         CC(1,1,-i)= b(7)
 947         CC(2,2,-i)=-b(7)
 948         CC(2,1,-i)=-b(9)
 949         CC(1,2,-i)=-b(9)
 950 c        CC(1,1,i)=0.0d0
 951 c        CC(2,2,i)=0.0d0
 952 c        CC(2,1,i)=0.0d0
 953 c        CC(1,2,i)=0.0d0
 954         Ctilde(1,1,i)=b(7)
 955         Ctilde(1,2,i)=b(9)
 956         Ctilde(2,1,i)=-b(9)
 957         Ctilde(2,2,i)=b(7)
 958         Ctilde(1,1,-i)=b(7)
 959         Ctilde(1,2,-i)=-b(9)
 960         Ctilde(2,1,-i)=b(9)
 961         Ctilde(2,2,-i)=b(7)
 962
 963 c        Ctilde(1,1,i)=0.0d0
 964 c        Ctilde(1,2,i)=0.0d0
 965 c        Ctilde(2,1,i)=0.0d0
 966 c        Ctilde(2,2,i)=0.0d0
 967         DD(1,1,i)= b(6)
 968         DD(2,2,i)=-b(6)
 969         DD(2,1,i)= b(8)
 970         DD(1,2,i)= b(8)
 971         DD(1,1,-i)= b(6)
 972         DD(2,2,-i)=-b(6)
 973         DD(2,1,-i)=-b(8)
 974         DD(1,2,-i)=-b(8)
 975 c        DD(1,1,i)=0.0d0
 976 c        DD(2,2,i)=0.0d0
 977 c        DD(2,1,i)=0.0d0
 978 c        DD(1,2,i)=0.0d0
 979         Dtilde(1,1,i)=b(6)
 980         Dtilde(1,2,i)=b(8)
 981         Dtilde(2,1,i)=-b(8)
 982         Dtilde(2,2,i)=b(6)
 983         Dtilde(1,1,-i)=b(6)
 984         Dtilde(1,2,-i)=-b(8)
 985         Dtilde(2,1,-i)=b(8)
 986         Dtilde(2,2,-i)=b(6)
 987
 988 c        Dtilde(1,1,i)=0.0d0
 989 c        Dtilde(1,2,i)=0.0d0
 990 c        Dtilde(2,1,i)=0.0d0
 991 c        Dtilde(2,2,i)=0.0d0
 992         EE(1,1,i)= b(10)+b(11)
 993         EE(2,2,i)=-b(10)+b(11)
 994         EE(2,1,i)= b(12)-b(13)
 995         EE(1,2,i)= b(12)+b(13)
 996         EE(1,1,-i)= b(10)+b(11)
 997         EE(2,2,-i)=-b(10)+b(11)
 998         EE(2,1,-i)=-b(12)+b(13)
 999         EE(1,2,-i)=-b(12)-b(13)
1000
1001 c        ee(1,1,i)=1.0d0
1002 c        ee(2,2,i)=1.0d0
1003 c        ee(2,1,i)=0.0d0
1004 c        ee(1,2,i)=0.0d0
1005 c        ee(2,1,i)=ee(1,2,i)
1006       enddo
1007       if (lprint) then
1008       do i=1,nloctyp
1009         write (iout,*) 'Type',i
1010         write (iout,*) 'B1'
1011         write(iout,*) B1(1,i),B1(2,i)
1012         write (iout,*) 'B2'
1013         write(iout,*) B2(1,i),B2(2,i)
1014         write (iout,*) 'CC'
1015         do j=1,2
1016           write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
1017         enddo
1018         write(iout,*) 'DD'
1019         do j=1,2
1020           write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
1021         enddo
1022         write(iout,*) 'EE'
1023         do j=1,2
1024           write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
1025         enddo
1026       enddo
1027       endif
1028
1029 C
1030 C Read electrostatic-interaction parameters
1031 C
1032       if (lprint) then
1033         write (iout,*)
1034         write (iout,'(/a)') 'Electrostatic interaction constants:'
1035         write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
1036      &            'IT','JT','APP','BPP','AEL6','AEL3'
1037       endif
1038       read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
1039       read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
1040       read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
1041       read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
1042       close (ielep)
1043       do i=1,2
1044         do j=1,2
1045         rri=rpp(i,j)**6
1046         app (i,j)=epp(i,j)*rri*rri
1047         bpp (i,j)=-2.0D0*epp(i,j)*rri
1048         ael6(i,j)=elpp6(i,j)*4.2D0**6
1049         ael3(i,j)=elpp3(i,j)*4.2D0**3
1050 c        lprint=.true.
1051         if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
1052      &                    ael6(i,j),ael3(i,j)
1053 c        lprint=.false.
1054         enddo
1055       enddo
1056 C
1057 C Read side-chain interaction parameters.
1058 C
1059       read (isidep,*,end=117,err=117) ipot,expon
1060       if (ipot.lt.1 .or. ipot.gt.5) then
1061         write (iout,'(2a)') 'Error while reading SC interaction',
1062      &               'potential file - unknown potential type.'
1063 #ifdef MPI
1064         call MPI_Finalize(Ierror)
1065 #endif
1066         stop
1067       endif
1068       expon2=expon/2
1069       if(me.eq.king)
1070      & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
1071      & ', exponents are ',expon,2*expon
1072       goto (10,20,30,30,40) ipot
1073 C----------------------- LJ potential ---------------------------------
1074    10 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
1075      &   (sigma0(i),i=1,ntyp)
1076       if (lprint) then
1077         write (iout,'(/a/)') 'Parameters of the LJ potential:'
1078         write (iout,'(a/)') 'The epsilon array:'
1079         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1080         write (iout,'(/a)') 'One-body parameters:'
1081         write (iout,'(a,4x,a)') 'residue','sigma'
1082         write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1083       endif
1084       goto 50
1085 C----------------------- LJK potential --------------------------------
1086    20 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
1087      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1088       if (lprint) then
1089         write (iout,'(/a/)') 'Parameters of the LJK potential:'
1090         write (iout,'(a/)') 'The epsilon array:'
1091         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1092         write (iout,'(/a)') 'One-body parameters:'
1093         write (iout,'(a,4x,2a)') 'residue','   sigma  ','    r0    '
1094         write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
1095      &        i=1,ntyp)
1096       endif
1097       goto 50
1098 C---------------------- GB or BP potential -----------------------------
1099    30 do i=1,ntyp
1100        read (isidep,*,end=116,err=116)(eps(i,j),j=i,ntyp)
1101       enddo
1102       read (isidep,*,end=116,err=116)(sigma0(i),i=1,ntyp)
1103       read (isidep,*,end=116,err=116)(sigii(i),i=1,ntyp)
1104       read (isidep,*,end=116,err=116)(chip(i),i=1,ntyp)
1105       read (isidep,*,end=116,err=116)(alp(i),i=1,ntyp)
1106 C now we start reading lipid
1107       do i=1,ntyp
1108        read (isidep,*,end=1161,err=1161)(epslip(i,j),j=i,ntyp)
1109        print *,"WARNING!!"
1110        do j=1,ntyp
1111        epslip(i,j)=epslip(i,j)+0.5d0
1112        enddo
1113       enddo
1114 C For the GB potential convert sigma'**2 into chi'
1115       if (ipot.eq.4) then
1116         do i=1,ntyp
1117           chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1118         enddo
1119       endif
1120       if (lprint) then
1121         write (iout,'(/a/)') 'Parameters of the BP potential:'
1122         write (iout,'(a/)') 'The epsilon array:'
1123         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1124         write (iout,'(/a)') 'One-body parameters:'
1125         write (iout,'(a,4x,4a)') 'residue','   sigma  ','s||/s_|_^2',
1126      &       '    chip  ','    alph  '
1127         write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1128      &                     chip(i),alp(i),i=1,ntyp)
1129       endif
1130       goto 50
1131 C--------------------- GBV potential -----------------------------------
1132    40 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
1133      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1134      &  (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1135       if (lprint) then
1136         write (iout,'(/a/)') 'Parameters of the GBV potential:'
1137         write (iout,'(a/)') 'The epsilon array:'
1138         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1139         write (iout,'(/a)') 'One-body parameters:'
1140         write (iout,'(a,4x,5a)') 'residue','   sigma  ','    r0    ',
1141      &      's||/s_|_^2','    chip  ','    alph  '
1142         write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1143      &           sigii(i),chip(i),alp(i),i=1,ntyp)
1144       endif
1145    50 continue
1146       close (isidep)
1147 C-----------------------------------------------------------------------
1148 C Calculate the "working" parameters of SC interactions.
1149       do i=2,ntyp
1150         do j=1,i-1
1151           eps(i,j)=eps(j,i)
1152         enddo
1153       enddo
1154       do i=1,ntyp
1155         do j=i,ntyp
1156           sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1157           sigma(j,i)=sigma(i,j)
1158           rs0(i,j)=dwa16*sigma(i,j)
1159           rs0(j,i)=rs0(i,j)
1160         enddo
1161       enddo
1162       if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1163      & 'Working parameters of the SC interactions:',
1164      & '     a    ','     b    ','   augm   ','  sigma ','   r0   ',
1165      & '  chi1   ','   chi2   '
1166       do i=1,ntyp
1167         do j=i,ntyp
1168           epsij=eps(i,j)
1169           if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1170             rrij=sigma(i,j)
1171           else
1172             rrij=rr0(i)+rr0(j)
1173           endif
1174           r0(i,j)=rrij
1175           r0(j,i)=rrij
1176           rrij=rrij**expon
1177           epsij=eps(i,j)
1178           sigeps=dsign(1.0D0,epsij)
1179           epsij=dabs(epsij)
1180           aa_aq(i,j)=epsij*rrij*rrij
1181           bb_aq(i,j)=-sigeps*epsij*rrij
1182           aa_aq(j,i)=aa_aq(i,j)
1183           bb_aq(j,i)=bb_aq(i,j)
1184           epsijlip=epslip(i,j)
1185           sigeps=dsign(1.0D0,epsijlip)
1186           epsijlip=dabs(epsijlip)
1187           aa_lip(i,j)=epsijlip*rrij*rrij
1188           bb_lip(i,j)=-sigeps*epsijlip*rrij
1189           aa_lip(j,i)=aa_lip(i,j)
1190           bb_lip(j,i)=bb_lip(i,j)
1191           if (ipot.gt.2) then
1192             sigt1sq=sigma0(i)**2
1193             sigt2sq=sigma0(j)**2
1194             sigii1=sigii(i)
1195             sigii2=sigii(j)
1196             ratsig1=sigt2sq/sigt1sq
1197             ratsig2=1.0D0/ratsig1
1198             chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1199             if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1200             rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1201           else
1202             rsum_max=sigma(i,j)
1203           endif
1204 c         if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1205             sigmaii(i,j)=rsum_max
1206             sigmaii(j,i)=rsum_max
1207 c         else
1208 c           sigmaii(i,j)=r0(i,j)
1209 c           sigmaii(j,i)=r0(i,j)
1210 c         endif
1211 cd        write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1212           if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1213             r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1214             augm(i,j)=epsij*r_augm**(2*expon)
1215 c           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1216             augm(j,i)=augm(i,j)
1217           else
1218             augm(i,j)=0.0D0
1219             augm(j,i)=0.0D0
1220           endif
1221           if (lprint) then
1222             write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1223      &      restyp(i),restyp(j),aa_aq(i,j),bb_aq(i,j),augm(i,j),
1224      &      sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1225           endif
1226         enddo
1227       enddo
1228 #ifdef OLDSCP
1229 C
1230 C Define the SC-p interaction constants (hard-coded; old style)
1231 C
1232       do i=1,ntyp
1233 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1234 C helix formation)
1235 c       aad(i,1)=0.3D0*4.0D0**12
1236 C Following line for constants currently implemented
1237 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1238         aad(i,1)=1.5D0*4.0D0**12
1239 c       aad(i,1)=0.17D0*5.6D0**12
1240         aad(i,2)=aad(i,1)
1241 C "Soft" SC-p repulsion
1242         bad(i,1)=0.0D0
1243 C Following line for constants currently implemented
1244 c       aad(i,1)=0.3D0*4.0D0**6
1245 C "Hard" SC-p repulsion
1246         bad(i,1)=3.0D0*4.0D0**6
1247 c       bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1248         bad(i,2)=bad(i,1)
1249 c       aad(i,1)=0.0D0
1250 c       aad(i,2)=0.0D0
1251 c       bad(i,1)=1228.8D0
1252 c       bad(i,2)=1228.8D0
1253       enddo
1254 #else
1255 C
1256 C 8/9/01 Read the SC-p interaction constants from file
1257 C
1258       do i=1,ntyp
1259         read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1260       enddo
1261       do i=1,ntyp
1262         aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1263         aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1264         bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1265         bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1266       enddo
1267 c      lprint=.true.
1268       if (lprint) then
1269         write (iout,*) "Parameters of SC-p interactions:"
1270         do i=1,ntyp
1271           write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1272      &     eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1273         enddo
1274       endif
1275 c      lprint=.false.
1276 #endif
1277 C
1278 C Define the constants of the disulfide bridge
1279 C
1280       ebr=-5.50D0
1281 c
1282 c Old arbitrary potential - commented out.
1283 c
1284 c      dbr= 4.20D0
1285 c      fbr= 3.30D0
1286 c
1287 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1288 c energy surface of diethyl disulfide.
1289 c A. Liwo and U. Kozlowska, 11/24/03
1290 c
1291       D0CM = 3.78d0
1292       AKCM = 15.1d0
1293       AKTH = 11.0d0
1294       AKCT = 12.0d0
1295       V1SS =-1.08d0
1296       V2SS = 7.61d0
1297       V3SS = 13.7d0
1298 c      akcm=0.0d0
1299 c      akth=0.0d0
1300 c      akct=0.0d0
1301 c      v1ss=0.0d0
1302 c      v2ss=0.0d0
1303 c      v3ss=0.0d0
1304
1305       if(me.eq.king) then
1306       write (iout,'(/a)') "Disulfide bridge parameters:"
1307       write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1308       write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1309       write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1310       write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1311      &  ' v3ss:',v3ss
1312       endif
1313       return
1314   111 write (iout,*) "Error reading bending energy parameters."
1315       goto 999
1316   112 write (iout,*) "Error reading rotamer energy parameters."
1317       goto 999
1318   113 write (iout,*) "Error reading torsional energy parameters."
1319       goto 999
1320   114 write (iout,*) "Error reading double torsional energy parameters."
1321       goto 999
1322   115 write (iout,*)
1323      &  "Error reading cumulant (multibody energy) parameters."
1324       goto 999
1325   116 write (iout,*) "Error reading electrostatic energy parameters."
1326       goto 999
1327  1161 write (iout,*) "Error reading electrostatic energy parameters.Lip"
1328       goto 999
1329   117 write (iout,*) "Error reading side chain interaction parameters."
1330       goto 999
1331   118 write (iout,*) "Error reading SCp interaction parameters."
1332       goto 999
1333   119 write (iout,*) "Error reading SCCOR parameters"
1334   999 continue
1335 #ifdef MPI
1336       call MPI_Finalize(Ierror)
1337 #endif
1338       stop
1339       return
1340       end
1341
1342
1343       subroutine getenv_loc(var, val)
1344       character(*) var, val
1345
1346 #ifdef WINIFL
1347       character(2000) line
1348       external ilen
1349
1350       open (196,file='env',status='old',readonly,shared)
1351       iread=0
1352 c      write(*,*)'looking for ',var
1353 10    read(196,*,err=11,end=11)line
1354       iread=index(line,var)
1355 c      write(*,*)iread,' ',var,' ',line
1356       if (iread.eq.0) go to 10
1357 c      write(*,*)'---> ',line
1358 11    continue
1359       if(iread.eq.0) then
1360 c       write(*,*)'CHUJ'
1361        val=''
1362       else
1363        iread=iread+ilen(var)+1
1364        read (line(iread:),*,err=12,end=12) val
1365 c       write(*,*)'OK: ',var,' = ',val
1366       endif
1367       close(196)
1368       return
1369 12    val=''
1370       close(196)
1371 #elif (defined CRAY)
1372       integer lennam,lenval,ierror
1373 c
1374 c        getenv using a POSIX call, useful on the T3D
1375 c        Sept 1996, comment out error check on advice of H. Pritchard
1376 c
1377       lennam = len(var)
1378       if(lennam.le.0) stop '--error calling getenv--'
1379       call pxfgetenv(var,lennam,val,lenval,ierror)
1380 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1381 #else
1382       call getenv(var,val)
1383 #endif
1384
1385       return
1386       end