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