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