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 For the GB potential convert sigma'**2 into chi'
1107       if (ipot.eq.4) then
1108         do i=1,ntyp
1109           chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1110         enddo
1111       endif
1112       if (lprint) then
1113         write (iout,'(/a/)') 'Parameters of the BP potential:'
1114         write (iout,'(a/)') 'The epsilon array:'
1115         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1116         write (iout,'(/a)') 'One-body parameters:'
1117         write (iout,'(a,4x,4a)') 'residue','   sigma  ','s||/s_|_^2',
1118      &       '    chip  ','    alph  '
1119         write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1120      &                     chip(i),alp(i),i=1,ntyp)
1121       endif
1122       goto 50
1123 C--------------------- GBV potential -----------------------------------
1124    40 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
1125      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1126      &  (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1127       if (lprint) then
1128         write (iout,'(/a/)') 'Parameters of the GBV potential:'
1129         write (iout,'(a/)') 'The epsilon array:'
1130         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1131         write (iout,'(/a)') 'One-body parameters:'
1132         write (iout,'(a,4x,5a)') 'residue','   sigma  ','    r0    ',
1133      &      's||/s_|_^2','    chip  ','    alph  '
1134         write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1135      &           sigii(i),chip(i),alp(i),i=1,ntyp)
1136       endif
1137    50 continue
1138       close (isidep)
1139 C-----------------------------------------------------------------------
1140 C Calculate the "working" parameters of SC interactions.
1141       do i=2,ntyp
1142         do j=1,i-1
1143           eps(i,j)=eps(j,i)
1144         enddo
1145       enddo
1146       do i=1,ntyp
1147         do j=i,ntyp
1148           sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1149           sigma(j,i)=sigma(i,j)
1150           rs0(i,j)=dwa16*sigma(i,j)
1151           rs0(j,i)=rs0(i,j)
1152         enddo
1153       enddo
1154       if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1155      & 'Working parameters of the SC interactions:',
1156      & '     a    ','     b    ','   augm   ','  sigma ','   r0   ',
1157      & '  chi1   ','   chi2   '
1158       do i=1,ntyp
1159         do j=i,ntyp
1160           epsij=eps(i,j)
1161           if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1162             rrij=sigma(i,j)
1163           else
1164             rrij=rr0(i)+rr0(j)
1165           endif
1166           r0(i,j)=rrij
1167           r0(j,i)=rrij
1168           rrij=rrij**expon
1169           epsij=eps(i,j)
1170           sigeps=dsign(1.0D0,epsij)
1171           epsij=dabs(epsij)
1172           aa(i,j)=epsij*rrij*rrij
1173           bb(i,j)=-sigeps*epsij*rrij
1174           aa(j,i)=aa(i,j)
1175           bb(j,i)=bb(i,j)
1176           if (ipot.gt.2) then
1177             sigt1sq=sigma0(i)**2
1178             sigt2sq=sigma0(j)**2
1179             sigii1=sigii(i)
1180             sigii2=sigii(j)
1181             ratsig1=sigt2sq/sigt1sq
1182             ratsig2=1.0D0/ratsig1
1183             chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1184             if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1185             rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1186           else
1187             rsum_max=sigma(i,j)
1188           endif
1189 c         if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1190             sigmaii(i,j)=rsum_max
1191             sigmaii(j,i)=rsum_max
1192 c         else
1193 c           sigmaii(i,j)=r0(i,j)
1194 c           sigmaii(j,i)=r0(i,j)
1195 c         endif
1196 cd        write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1197           if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1198             r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1199             augm(i,j)=epsij*r_augm**(2*expon)
1200 c           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1201             augm(j,i)=augm(i,j)
1202           else
1203             augm(i,j)=0.0D0
1204             augm(j,i)=0.0D0
1205           endif
1206           if (lprint) then
1207             write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1208      &      restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1209      &      sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1210           endif
1211         enddo
1212       enddo
1213 #ifdef OLDSCP
1214 C
1215 C Define the SC-p interaction constants (hard-coded; old style)
1216 C
1217       do i=1,ntyp
1218 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1219 C helix formation)
1220 c       aad(i,1)=0.3D0*4.0D0**12
1221 C Following line for constants currently implemented
1222 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1223         aad(i,1)=1.5D0*4.0D0**12
1224 c       aad(i,1)=0.17D0*5.6D0**12
1225         aad(i,2)=aad(i,1)
1226 C "Soft" SC-p repulsion
1227         bad(i,1)=0.0D0
1228 C Following line for constants currently implemented
1229 c       aad(i,1)=0.3D0*4.0D0**6
1230 C "Hard" SC-p repulsion
1231         bad(i,1)=3.0D0*4.0D0**6
1232 c       bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1233         bad(i,2)=bad(i,1)
1234 c       aad(i,1)=0.0D0
1235 c       aad(i,2)=0.0D0
1236 c       bad(i,1)=1228.8D0
1237 c       bad(i,2)=1228.8D0
1238       enddo
1239 #else
1240 C
1241 C 8/9/01 Read the SC-p interaction constants from file
1242 C
1243       do i=1,ntyp
1244         read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1245       enddo
1246       do i=1,ntyp
1247         aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1248         aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1249         bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1250         bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1251       enddo
1252 c      lprint=.true.
1253       if (lprint) then
1254         write (iout,*) "Parameters of SC-p interactions:"
1255         do i=1,ntyp
1256           write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1257      &     eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1258         enddo
1259       endif
1260 c      lprint=.false.
1261 #endif
1262 C
1263 C Define the constants of the disulfide bridge
1264 C
1265       ebr=-5.50D0
1266 c
1267 c Old arbitrary potential - commented out.
1268 c
1269 c      dbr= 4.20D0
1270 c      fbr= 3.30D0
1271 c
1272 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1273 c energy surface of diethyl disulfide.
1274 c A. Liwo and U. Kozlowska, 11/24/03
1275 c
1276       D0CM = 3.78d0
1277       AKCM = 15.1d0
1278       AKTH = 11.0d0
1279       AKCT = 12.0d0
1280       V1SS =-1.08d0
1281       V2SS = 7.61d0
1282       V3SS = 13.7d0
1283 c      akcm=0.0d0
1284 c      akth=0.0d0
1285 c      akct=0.0d0
1286 c      v1ss=0.0d0
1287 c      v2ss=0.0d0
1288 c      v3ss=0.0d0
1289
1290       if(me.eq.king) then
1291       write (iout,'(/a)') "Disulfide bridge parameters:"
1292       write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1293       write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1294       write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1295       write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1296      &  ' v3ss:',v3ss
1297       endif
1298       return
1299   111 write (iout,*) "Error reading bending energy parameters."
1300       goto 999
1301   112 write (iout,*) "Error reading rotamer energy parameters."
1302       goto 999
1303   113 write (iout,*) "Error reading torsional energy parameters."
1304       goto 999
1305   114 write (iout,*) "Error reading double torsional energy parameters."
1306       goto 999
1307   115 write (iout,*)
1308      &  "Error reading cumulant (multibody energy) parameters."
1309       goto 999
1310   116 write (iout,*) "Error reading electrostatic energy parameters."
1311       goto 999
1312   117 write (iout,*) "Error reading side chain interaction parameters."
1313       goto 999
1314   118 write (iout,*) "Error reading SCp interaction parameters."
1315       goto 999
1316   119 write (iout,*) "Error reading SCCOR parameters"
1317   999 continue
1318 #ifdef MPI
1319       call MPI_Finalize(Ierror)
1320 #endif
1321       stop
1322       return
1323       end
1324
1325
1326       subroutine getenv_loc(var, val)
1327       character(*) var, val
1328
1329 #ifdef WINIFL
1330       character(2000) line
1331       external ilen
1332
1333       open (196,file='env',status='old',readonly,shared)
1334       iread=0
1335 c      write(*,*)'looking for ',var
1336 10    read(196,*,err=11,end=11)line
1337       iread=index(line,var)
1338 c      write(*,*)iread,' ',var,' ',line
1339       if (iread.eq.0) go to 10
1340 c      write(*,*)'---> ',line
1341 11    continue
1342       if(iread.eq.0) then
1343 c       write(*,*)'CHUJ'
1344        val=''
1345       else
1346        iread=iread+ilen(var)+1
1347        read (line(iread:),*,err=12,end=12) val
1348 c       write(*,*)'OK: ',var,' = ',val
1349       endif
1350       close(196)
1351       return
1352 12    val=''
1353       close(196)
1354 #elif (defined CRAY)
1355       integer lennam,lenval,ierror
1356 c
1357 c        getenv using a POSIX call, useful on the T3D
1358 c        Sept 1996, comment out error check on advice of H. Pritchard
1359 c
1360       lennam = len(var)
1361       if(lennam.le.0) stop '--error calling getenv--'
1362       call pxfgetenv(var,lennam,val,lenval,ierror)
1363 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1364 #else
1365       call getenv(var,val)
1366 #endif
1367
1368       return
1369       end