source/wham/src-M/parmread.F

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