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