source/wham/src-M-homology/parmread.F

   1       subroutine parmread(iparm,*)
   2 C
   3 C Read the parameters of the probability distributions of the virtual-bond
   4 C valence angles and the side chains and energy parameters.
   5 C
   6       implicit real*8 (a-h,o-z)
   7       include 'DIMENSIONS'
   8       include 'DIMENSIONS.ZSCOPT'
   9       include 'DIMENSIONS.FREE'
  10       include 'COMMON.IOUNITS'
  11       include 'COMMON.CHAIN'
  12       include 'COMMON.INTERACT'
  13       include 'COMMON.GEO'
  14       include 'COMMON.LOCAL'
  15       include 'COMMON.TORSION'
  16       include 'COMMON.FFIELD'
  17       include 'COMMON.NAMES'
  18       include 'COMMON.SBRIDGE'
  19       include 'COMMON.WEIGHTS'
  20       include 'COMMON.ENEPS'
  21       include 'COMMON.SCCOR'
  22       include 'COMMON.SCROT'
  23       include 'COMMON.FREE'
  24       include 'COMMON.CONTROL'
  25       include 'COMMON.LANGEVIN'
  26       character*1 t1,t2,t3
  27       character*1 onelett(4) /"G","A","P","D"/
  28       character*1 toronelet(-2:2) /"p","a","G","A","P"/
  29       logical lprint
  30       dimension blower(3,3,maxlob)
  31       character*800 controlcard
  32       character*256 bondname_t,thetname_t,rotname_t,torname_t,
  33      & tordname_t,fouriername_t,elename_t,sidename_t,scpname_t,
  34      & sccorname_t
  35       integer ilen
  36       external ilen
  37       character*16 key
  38       integer iparm
  39       character*6 res1
  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       wsaxs=ww(26)
 144
 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,mp,ip,pstok
 217       do i=1,ntyp
 218         nbondterm(i)=1
 219         read (ibond,*) vbldsc0(1,i),aksc(1,i),msc(i),isc(i),restok(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,mp,ip,pstok
 229       do i=1,ntyp
 230         read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),
 231      &   j=1,nbondterm(i)),msc(i),isc(i),restok(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       read (ithep,*) nthetyp,ntheterm,ntheterm2,
 362      &  ntheterm3,nsingle,ndouble
 363       nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
 364       read (ithep,*) (ithetyp(i),i=1,ntyp1)
 365       do i=-ntyp1,-1
 366         ithetyp(i)=-ithetyp(-i)
 367       enddo
 368 c      write (iout,*) "tu dochodze"
 369       do iblock=1,2
 370       do i=-maxthetyp,maxthetyp
 371         do j=-maxthetyp,maxthetyp
 372           do k=-maxthetyp,maxthetyp
 373             aa0thet(i,j,k,iblock)=0.0d0
 374             do l=1,ntheterm
 375               aathet(l,i,j,k,iblock)=0.0d0
 376             enddo
 377             do l=1,ntheterm2
 378               do m=1,nsingle
 379                 bbthet(m,l,i,j,k,iblock)=0.0d0
 380                 ccthet(m,l,i,j,k,iblock)=0.0d0
 381                 ddthet(m,l,i,j,k,iblock)=0.0d0
 382                 eethet(m,l,i,j,k,iblock)=0.0d0
 383               enddo
 384             enddo
 385             do l=1,ntheterm3
 386               do m=1,ndouble
 387                 do mm=1,ndouble
 388                  ffthet(mm,m,l,i,j,k,iblock)=0.0d0
 389                  ggthet(mm,m,l,i,j,k,iblock)=0.0d0
 390                 enddo
 391               enddo
 392             enddo
 393           enddo
 394         enddo
 395       enddo
 396       enddo
 397 C      write (iout,*) "KURWA1"
 398       do iblock=1,2
 399       do i=0,nthetyp
 400         do j=-nthetyp,nthetyp
 401           do k=-nthetyp,nthetyp
 402             read (ithep,'(6a)') res1
 403 c            write(iout,*) res1,i,j,k
 404             read (ithep,*) aa0thet(i,j,k,iblock)
 405             read (ithep,*)(aathet(l,i,j,k,iblock),l=1,ntheterm)
 406             read (ithep,*)
 407      &       ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
 408      &        (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
 409      &        (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
 410      &        (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle)
 411      &        ,ll=1,ntheterm2)
 412             read (ithep,*)
 413      &      (((ffthet(llll,lll,ll,i,j,k,iblock),
 414      &      ffthet(lll,llll,ll,i,j,k,iblock),
 415      &         ggthet(llll,lll,ll,i,j,k,iblock)
 416      &        ,ggthet(lll,llll,ll,i,j,k,iblock),
 417      &         llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
 418           enddo
 419         enddo
 420       enddo
 421 C
 422 C For dummy ends assign glycine-type coefficients of theta-only terms; the
 423 C coefficients of theta-and-gamma-dependent terms are zero.
 424 C
 425       do i=1,nthetyp
 426         do j=1,nthetyp
 427           do l=1,ntheterm
 428             aathet(l,i,j,nthetyp+1,iblock)=0.0d0
 429             aathet(l,nthetyp+1,i,j,iblock)=0.0d0
 430           enddo
 431           aa0thet(i,j,nthetyp+1,iblock)=0.0d0
 432           aa0thet(nthetyp+1,i,j,iblock)=0.0d0
 433         enddo
 434         do l=1,ntheterm
 435           aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
 436         enddo
 437         aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
 438       enddo
 439       enddo
 440 C       write(iout,*) "KURWA1.5"
 441 C Substitution for D aminoacids from symmetry.
 442       do iblock=1,2
 443       do i=-nthetyp,0
 444         do j=-nthetyp,nthetyp
 445           do k=-nthetyp,nthetyp
 446            aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
 447            do l=1,ntheterm
 448            aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
 449            enddo
 450            do ll=1,ntheterm2
 451             do lll=1,nsingle
 452             bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
 453             ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
 454             ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
 455             eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
 456             enddo
 457           enddo
 458           do ll=1,ntheterm3
 459            do lll=2,ndouble
 460             do llll=1,lll-1
 461             ffthet(llll,lll,ll,i,j,k,iblock)=
 462      &      ffthet(llll,lll,ll,-i,-j,-k,iblock)
 463             ffthet(lll,llll,ll,i,j,k,iblock)=
 464      &      ffthet(lll,llll,ll,-i,-j,-k,iblock)
 465             ggthet(llll,lll,ll,i,j,k,iblock)=
 466      &      -ggthet(llll,lll,ll,-i,-j,-k,iblock)
 467             ggthet(lll,llll,ll,i,j,k,iblock)=
 468      &      -ggthet(lll,llll,ll,-i,-j,-k,iblock)
 469             enddo !ll
 470            enddo  !lll
 471           enddo   !llll
 472          enddo    !k
 473         enddo     !j
 474        enddo      !i
 475       enddo       !iblock
 476
 477 C
 478 C Control printout of the coefficients of virtual-bond-angle potentials
 479 C
 480       if (lprint) then
 481         write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
 482         do i=1,nthetyp+1
 483           do j=1,nthetyp+1
 484             do k=1,nthetyp+1
 485               write (iout,'(//4a)')
 486      &         'Type ',onelett(i),onelett(j),onelett(k)
 487               write (iout,'(//a,10x,a)') " l","a[l]"
 488               write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
 489               write (iout,'(i2,1pe15.5)')
 490      &           (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
 491             do l=1,ntheterm2
 492               write (iout,'(//2h m,4(9x,a,3h[m,i1,1h]))')
 493      &          "b",l,"c",l,"d",l,"e",l
 494               do m=1,nsingle
 495                 write (iout,'(i2,4(1pe15.5))') m,
 496      &          bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),
 497      &          ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
 498               enddo
 499             enddo
 500             do l=1,ntheterm3
 501               write (iout,'(//3hm,n,4(6x,a,5h[m,n,i1,1h]))')
 502      &          "f+",l,"f-",l,"g+",l,"g-",l
 503               do m=2,ndouble
 504                 do n=1,m-1
 505                   write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
 506      &              ffthet(n,m,l,i,j,k,iblock),
 507      &              ffthet(m,n,l,i,j,k,iblock),
 508      &              ggthet(n,m,l,i,j,k,iblock),
 509      &              ggthet(m,n,l,i,j,k,iblock)
 510                 enddo
 511               enddo
 512             enddo
 513           enddo
 514         enddo
 515       enddo
 516       call flush(iout)
 517       endif
 518 #endif
 519
 520 #ifdef CRYST_SC
 521 C
 522 C Read the parameters of the probability distribution/energy expression
 523 C of the side chains.
 524 C
 525       do i=1,ntyp
 526         read (irotam,'(3x,i3,f8.3)') nlob(i),dsc(i)
 527         if (i.eq.10) then
 528           dsc_inv(i)=0.0D0
 529         else
 530           dsc_inv(i)=1.0D0/dsc(i)
 531         endif
 532         if (i.ne.10) then
 533         do j=1,nlob(i)
 534           do k=1,3
 535             do l=1,3
 536               blower(l,k,j)=0.0D0
 537             enddo
 538           enddo
 539         enddo
 540         bsc(1,i)=0.0D0
 541         read(irotam,*)(censc(k,1,i),k=1,3),((blower(k,l,1),l=1,k),k=1,3)
 542         censc(1,1,-i)=censc(1,1,i)
 543         censc(2,1,-i)=censc(2,1,i)
 544         censc(3,1,-i)=-censc(3,1,i)
 545         do j=2,nlob(i)
 546           read (irotam,*) bsc(j,i)
 547           read (irotam,*) (censc(k,j,i),k=1,3),
 548      &                                 ((blower(k,l,j),l=1,k),k=1,3)
 549         censc(1,j,-i)=censc(1,j,i)
 550         censc(2,j,-i)=censc(2,j,i)
 551         censc(3,j,-i)=-censc(3,j,i)
 552 C BSC is amplitude of Gaussian
 553         enddo
 554         do j=1,nlob(i)
 555           do k=1,3
 556             do l=1,k
 557               akl=0.0D0
 558               do m=1,3
 559                 akl=akl+blower(k,m,j)*blower(l,m,j)
 560               enddo
 561               gaussc(k,l,j,i)=akl
 562               gaussc(l,k,j,i)=akl
 563              if (((k.eq.3).and.(l.ne.3))
 564      &        .or.((l.eq.3).and.(k.ne.3))) then
 565                 gaussc(k,l,j,-i)=-akl
 566                 gaussc(l,k,j,-i)=-akl
 567               else
 568                 gaussc(k,l,j,-i)=akl
 569                 gaussc(l,k,j,-i)=akl
 570               endif
 571             enddo
 572           enddo
 573         enddo
 574         endif
 575       enddo
 576       close (irotam)
 577       if (lprint) then
 578         write (iout,'(/a)') 'Parameters of side-chain local geometry'
 579         do i=1,ntyp
 580           nlobi=nlob(i)
 581           if (nlobi.gt.0) then
 582           write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
 583      &     ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
 584 c          write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
 585 c          write (iout,'(a,f10.4,4(16x,f10.4))')
 586 c     &                             'Center  ',(bsc(j,i),j=1,nlobi)
 587 c          write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),j=1,nlobi)
 588            write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
 589      &                             'log h',(bsc(j,i),j=1,nlobi)
 590            write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
 591      &    'x',((censc(k,j,i),k=1,3),j=1,nlobi)
 592 c          write (iout,'(a)')
 593 c         do j=1,nlobi
 594 c           ind=0
 595 c           do k=1,3
 596 c             do l=1,k
 597 c              ind=ind+1
 598 c              blower(k,l,j)=gaussc(ind,j,i)
 599 c             enddo
 600 c           enddo
 601 c         enddo
 602           do k=1,3
 603             write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
 604      &                 ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
 605           enddo
 606           endif
 607         enddo
 608       endif
 609 #else
 610 C
 611 C Read scrot parameters for potentials determined from all-atom AM1 calculations
 612 C added by Urszula Kozlowska 07/11/2007
 613 C
 614       do i=1,ntyp
 615         read (irotam,*)
 616        if (i.eq.10) then
 617          read (irotam,*)
 618        else
 619          do j=1,65
 620            read(irotam,*) sc_parmin(j,i)
 621          enddo
 622        endif
 623       enddo
 624 #endif
 625       close(irotam)
 626 #ifdef CRYST_TOR
 627 C
 628 C Read torsional parameters in old format
 629 C
 630       read (itorp,*) ntortyp,nterm_old
 631       write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
 632       read (itorp,*) (itortyp(i),i=1,ntyp)
 633       do i=1,ntortyp
 634         do j=1,ntortyp
 635           read (itorp,'(a)')
 636           do k=1,nterm_old
 637             read (itorp,*) kk,v1(k,j,i),v2(k,j,i)
 638           enddo
 639         enddo
 640       enddo
 641       close (itorp)
 642       if (lprint) then
 643         write (iout,'(/a/)') 'Torsional constants:'
 644         do i=1,ntortyp
 645           do j=1,ntortyp
 646             write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
 647             write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
 648           enddo
 649         enddo
 650       endif
 651
 652
 653 #else
 654 C
 655 C Read torsional parameters
 656 C
 657       read (itorp,*) ntortyp
 658       read (itorp,*) (itortyp(i),i=1,ntyp)
 659       write (iout,*) 'ntortyp',ntortyp
 660       do iblock=1,2
 661       do i=-ntyp,-1
 662        itortyp(i)=-itortyp(-i)
 663       enddo
 664 c      write (iout,*) 'ntortyp',ntortyp
 665       do i=0,ntortyp-1
 666         do j=-ntortyp+1,ntortyp-1
 667           read (itorp,*) nterm(i,j,iblock),
 668      &          nlor(i,j,iblock)
 669           nterm(-i,-j,iblock)=nterm(i,j,iblock)
 670           nlor(-i,-j,iblock)=nlor(i,j,iblock)
 671           v0ij=0.0d0
 672           si=-1.0d0
 673           do k=1,nterm(i,j,iblock)
 674             read (itorp,*) kk,v1(k,i,j,iblock),
 675      &      v2(k,i,j,iblock)
 676             v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
 677             v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
 678             v0ij=v0ij+si*v1(k,i,j,iblock)
 679             si=-si
 680          enddo
 681           do k=1,nlor(i,j,iblock)
 682             read (itorp,*) kk,vlor1(k,i,j),
 683      &        vlor2(k,i,j),vlor3(k,i,j)
 684             v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
 685           enddo
 686           v0(i,j,iblock)=v0ij
 687           v0(-i,-j,iblock)=v0ij
 688         enddo
 689       enddo
 690       enddo
 691       close (itorp)
 692       if (lprint) then
 693         write (iout,'(/a/)') 'Torsional constants:'
 694         do i=1,ntortyp
 695           do j=1,ntortyp
 696             write (iout,*) 'ityp',i,' jtyp',j
 697             write (iout,*) 'Fourier constants'
 698             do k=1,nterm(i,j,iblock)
 699               write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
 700      &        v2(k,i,j,iblock)
 701             enddo
 702             write (iout,*) 'Lorenz constants'
 703             do k=1,nlor(i,j,iblock)
 704               write (iout,'(3(1pe15.5))')
 705      &         vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
 706             enddo
 707           enddo
 708         enddo
 709       endif
 710 C
 711 C 6/23/01 Read parameters for double torsionals
 712 C
 713       do iblock=1,2
 714       do i=0,ntortyp-1
 715         do j=-ntortyp+1,ntortyp-1
 716           do k=-ntortyp+1,ntortyp-1
 717             read (itordp,'(3a1)') t1,t2,t3
 718 c              write (iout,*) "OK onelett",
 719 c     &         i,j,k,t1,t2,t3
 720
 721             if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
 722      &        .or. t3.ne.toronelet(k)) then
 723               write (iout,*) "Error in double torsional parameter file",
 724      &         i,j,k,t1,t2,t3
 725 #ifdef MPI
 726               call MPI_Finalize(Ierror)
 727 #endif
 728                stop "Error in double torsional parameter file"
 729             endif
 730           read (itordp,*) ntermd_1(i,j,k,iblock),
 731      &         ntermd_2(i,j,k,iblock)
 732             ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
 733             ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
 734             read (itordp,*) (v1c(1,l,i,j,k,iblock),l=1,
 735      &         ntermd_1(i,j,k,iblock))
 736             read (itordp,*) (v1s(1,l,i,j,k,iblock),l=1,
 737      &         ntermd_1(i,j,k,iblock))
 738             read (itordp,*) (v1c(2,l,i,j,k,iblock),l=1,
 739      &         ntermd_1(i,j,k,iblock))
 740             read (itordp,*) (v1s(2,l,i,j,k,iblock),l=1,
 741      &         ntermd_1(i,j,k,iblock))
 742 C Martix of D parameters for one dimesional foureir series
 743             do l=1,ntermd_1(i,j,k,iblock)
 744              v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
 745              v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
 746              v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
 747              v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
 748 c            write(iout,*) "whcodze" ,
 749 c     & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
 750             enddo
 751             read (itordp,*) ((v2c(l,m,i,j,k,iblock),
 752      &         v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
 753      &         v2s(m,l,i,j,k,iblock),
 754      &         m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
 755 C Martix of D parameters for two dimesional fourier series
 756             do l=1,ntermd_2(i,j,k,iblock)
 757              do m=1,l-1
 758              v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
 759              v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
 760              v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
 761              v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
 762              enddo!m
 763             enddo!l
 764           enddo!k
 765         enddo!j
 766       enddo!i
 767       enddo!iblock
 768       if (lprint) then
 769       write (iout,*)
 770       write (iout,*) 'Constants for double torsionals'
 771       do iblock=1,2
 772       do i=0,ntortyp-1
 773         do j=-ntortyp+1,ntortyp-1
 774           do k=-ntortyp+1,ntortyp-1
 775             write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
 776      &        ' nsingle',ntermd_1(i,j,k,iblock),
 777      &        ' ndouble',ntermd_2(i,j,k,iblock)
 778             write (iout,*)
 779             write (iout,*) 'Single angles:'
 780             do l=1,ntermd_1(i,j,k,iblock)
 781               write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
 782      &           v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
 783      &           v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),
 784      &           v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
 785             enddo
 786             write (iout,*)
 787             write (iout,*) 'Pairs of angles:'
 788             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
 789             do l=1,ntermd_2(i,j,k,iblock)
 790               write (iout,'(i5,20f10.5)')
 791      &         l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
 792             enddo
 793             write (iout,*)
 794            write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
 795             do l=1,ntermd_2(i,j,k,iblock)
 796               write (iout,'(i5,20f10.5)')
 797      &         l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),
 798      &         (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
 799             enddo
 800             write (iout,*)
 801           enddo
 802         enddo
 803       enddo
 804       enddo
 805       endif
 806 #endif
 807 C Read of Side-chain backbone correlation parameters
 808 C Modified 11 May 2012 by Adasko
 809 CCC
 810 C
 811       read (isccor,*) nsccortyp
 812       read (isccor,*) (isccortyp(i),i=1,ntyp)
 813       do i=-ntyp,-1
 814         isccortyp(i)=-isccortyp(-i)
 815       enddo
 816       iscprol=isccortyp(20)
 817 c      write (iout,*) 'ntortyp',ntortyp
 818       maxinter=3
 819 cc maxinter is maximum interaction sites
 820       do l=1,maxinter
 821       do i=1,nsccortyp
 822         do j=1,nsccortyp
 823           read (isccor,*)
 824      &nterm_sccor(i,j),nlor_sccor(i,j)
 825 c          write (iout,*) nterm_sccor(i,j)
 826           v0ijsccor=0.0d0
 827           v0ijsccor1=0.0d0
 828           v0ijsccor2=0.0d0
 829           v0ijsccor3=0.0d0
 830           si=-1.0d0
 831           nterm_sccor(-i,j)=nterm_sccor(i,j)
 832           nterm_sccor(-i,-j)=nterm_sccor(i,j)
 833           nterm_sccor(i,-j)=nterm_sccor(i,j)
 834 c          write (iout,*) nterm_sccor(i,j),nterm_sccor(-i,j),
 835 c     &    nterm_sccor(-i,-j),nterm_sccor(i,-j)
 836           do k=1,nterm_sccor(i,j)
 837             read (isccor,*) kk,v1sccor(k,l,i,j)
 838      &    ,v2sccor(k,l,i,j)
 839             if (j.eq.iscprol) then
 840              if (i.eq.isccortyp(10)) then
 841              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
 842              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
 843              else
 844              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
 845      &                        +v2sccor(k,l,i,j)*dsqrt(0.75d0)
 846              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
 847      &                        +v1sccor(k,l,i,j)*dsqrt(0.75d0)
 848              v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
 849              v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
 850              v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
 851              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
 852              endif
 853             else
 854              if (i.eq.isccortyp(10)) then
 855              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
 856              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
 857              else
 858                if (j.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              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              v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
 867              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
 868                 endif
 869                endif
 870             endif
 871             v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
 872             v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
 873             v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
 874             v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
 875             si=-si
 876            enddo
 877           do k=1,nlor_sccor(i,j)
 878             read (isccor,*) kk,vlor1sccor(k,i,j),
 879      &        vlor2sccor(k,i,j),vlor3sccor(k,i,j)
 880             v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
 881      &(1+vlor3sccor(k,i,j)**2)
 882           enddo
 883           v0sccor(l,i,j)=v0ijsccor
 884           v0sccor(l,-i,j)=v0ijsccor1
 885           v0sccor(l,i,-j)=v0ijsccor2
 886           v0sccor(l,-i,-j)=v0ijsccor3
 887           enddo
 888         enddo
 889       enddo
 890       close (isccor)
 891       if (lprint) then
 892         write (iout,'(/a/)') 'Torsional constants of SCCORR:'
 893         do i=1,nsccortyp
 894           do j=1,nsccortyp
 895             write (iout,*) 'ityp',i,' jtyp',j
 896             write (iout,*) 'Fourier constants'
 897             do k=1,nterm_sccor(i,j)
 898               write (iout,'(2(1pe15.5))')
 899      & v1sccor(k,l,i,j),v2sccor(k,l,i,j)
 900             enddo
 901             write (iout,*) 'Lorenz constants'
 902             do k=1,nlor_sccor(i,j)
 903               write (iout,'(3(1pe15.5))')
 904      &         vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
 905             enddo
 906           enddo
 907         enddo
 908       endif
 909 C
 910 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
 911 C         interaction energy of the Gly, Ala, and Pro prototypes.
 912 C
 913       read (ifourier,*) nloctyp
 914       do i=0,nloctyp-1
 915         read (ifourier,*)
 916         read (ifourier,*) (b(ii,i),ii=1,13)
 917         if (lprint) then
 918         write (iout,*) 'Type',i
 919         write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
 920         endif
 921         B1(1,i)  = b(3,i)
 922         B1(2,i)  = b(5,i)
 923         B1(1,-i) = b(3,i)
 924         B1(2,-i) = -b(5,i)
 925 c        b1(1,i)=0.0d0
 926 c        b1(2,i)=0.0d0
 927         B1tilde(1,i) = b(3,i)
 928         B1tilde(2,i) =-b(5,i)
 929         B1tilde(1,-i) =-b(3,i)
 930         B1tilde(2,-i) =b(5,i)
 931 c        b1tilde(1,i)=0.0d0
 932 c        b1tilde(2,i)=0.0d0
 933         B2(1,i)  = b(2,i)
 934         B2(2,i)  = b(4,i)
 935         B2(1,-i)  =b(2,i)
 936         B2(2,-i)  =-b(4,i)
 937
 938 c        b2(1,i)=0.0d0
 939 c        b2(2,i)=0.0d0
 940         CC(1,1,i)= b(7,i)
 941         CC(2,2,i)=-b(7,i)
 942         CC(2,1,i)= b(9,i)
 943         CC(1,2,i)= b(9,i)
 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 c        CC(1,1,i)=0.0d0
 949 c        CC(2,2,i)=0.0d0
 950 c        CC(2,1,i)=0.0d0
 951 c        CC(1,2,i)=0.0d0
 952         Ctilde(1,1,i)=b(7,i)
 953         Ctilde(1,2,i)=b(9,i)
 954         Ctilde(2,1,i)=-b(9,i)
 955         Ctilde(2,2,i)=b(7,i)
 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
 961 c        Ctilde(1,1,i)=0.0d0
 962 c        Ctilde(1,2,i)=0.0d0
 963 c        Ctilde(2,1,i)=0.0d0
 964 c        Ctilde(2,2,i)=0.0d0
 965         DD(1,1,i)= b(6,i)
 966         DD(2,2,i)=-b(6,i)
 967         DD(2,1,i)= b(8,i)
 968         DD(1,2,i)= b(8,i)
 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 c        DD(1,1,i)=0.0d0
 974 c        DD(2,2,i)=0.0d0
 975 c        DD(2,1,i)=0.0d0
 976 c        DD(1,2,i)=0.0d0
 977         Dtilde(1,1,i)=b(6,i)
 978         Dtilde(1,2,i)=b(8,i)
 979         Dtilde(2,1,i)=-b(8,i)
 980         Dtilde(2,2,i)=b(6,i)
 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
 986 c        Dtilde(1,1,i)=0.0d0
 987 c        Dtilde(1,2,i)=0.0d0
 988 c        Dtilde(2,1,i)=0.0d0
 989 c        Dtilde(2,2,i)=0.0d0
 990         EE(1,1,i)= b(10,i)+b(11,i)
 991         EE(2,2,i)=-b(10,i)+b(11,i)
 992         EE(2,1,i)= b(12,i)-b(13,i)
 993         EE(1,2,i)= b(12,i)+b(13,i)
 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
 999 c        ee(1,1,i)=1.0d0
1000 c        ee(2,2,i)=1.0d0
1001 c        ee(2,1,i)=0.0d0
1002 c        ee(1,2,i)=0.0d0
1003 c        ee(2,1,i)=ee(1,2,i)
1004
1005       enddo
1006       if (lprint) then
1007       do i=1,nloctyp
1008         write (iout,*) 'Type',i
1009         write (iout,*) 'B1'
1010 c        write (iout,'(f10.5)') B1(:,i)
1011         write(iout,*) B1(1,i),B1(2,i)
1012         write (iout,*) 'B2'
1013 c        write (iout,'(f10.5)') B2(:,i)
1014         write(iout,*) B2(1,i),B2(2,i)
1015         write (iout,*) 'CC'
1016         do j=1,2
1017           write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
1018         enddo
1019         write(iout,*) 'DD'
1020         do j=1,2
1021           write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
1022         enddo
1023         write(iout,*) 'EE'
1024         do j=1,2
1025           write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
1026         enddo
1027       enddo
1028       endif
1029 C
1030 C Read electrostatic-interaction parameters
1031 C
1032       if (lprint) then
1033         write (iout,'(/a)') 'Electrostatic interaction constants:'
1034         write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
1035      &            'IT','JT','APP','BPP','AEL6','AEL3'
1036       endif
1037       read (ielep,*) ((epp(i,j),j=1,2),i=1,2)
1038       read (ielep,*) ((rpp(i,j),j=1,2),i=1,2)
1039       read (ielep,*) ((elpp6(i,j),j=1,2),i=1,2)
1040       read (ielep,*) ((elpp3(i,j),j=1,2),i=1,2)
1041       close (ielep)
1042       do i=1,2
1043         do j=1,2
1044         rri=rpp(i,j)**6
1045         app (i,j)=epp(i,j)*rri*rri
1046         bpp (i,j)=-2.0D0*epp(i,j)*rri
1047         ael6(i,j)=elpp6(i,j)*4.2D0**6
1048         ael3(i,j)=elpp3(i,j)*4.2D0**3
1049         lprint=.true.
1050         if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
1051      &                    ael6(i,j),ael3(i,j)
1052          lprint=.false.
1053         enddo
1054       enddo
1055 C
1056 C Read side-chain interaction parameters.
1057 C
1058       read (isidep,*) ipot,expon
1059       if (ipot.lt.1 .or. ipot.gt.5) then
1060         write (iout,'(2a)') 'Error while reading SC interaction',
1061      &               'potential file - unknown potential type.'
1062         stop
1063       endif
1064       expon2=expon/2
1065       write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
1066      & ', exponents are ',expon,2*expon
1067       goto (10,20,30,30,40) ipot
1068 C----------------------- LJ potential ---------------------------------
1069    10 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
1070       if (lprint) then
1071         write (iout,'(/a/)') 'Parameters of the LJ potential:'
1072         write (iout,'(a/)') 'The epsilon array:'
1073         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1074         write (iout,'(/a)') 'One-body parameters:'
1075         write (iout,'(a,4x,a)') 'residue','sigma'
1076         write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1077       endif
1078       goto 50
1079 C----------------------- LJK potential --------------------------------
1080    20 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
1081      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1082       if (lprint) then
1083         write (iout,'(/a/)') 'Parameters of the LJK potential:'
1084         write (iout,'(a/)') 'The epsilon array:'
1085         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1086         write (iout,'(/a)') 'One-body parameters:'
1087         write (iout,'(a,4x,2a)') 'residue','   sigma  ','    r0    '
1088         write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
1089      &        i=1,ntyp)
1090       endif
1091       goto 50
1092 C---------------------- GB or BP potential -----------------------------
1093    30 do i=1,ntyp
1094        read (isidep,*)(eps(i,j),j=i,ntyp)
1095       enddo
1096       read (isidep,*)(sigma0(i),i=1,ntyp)
1097       read (isidep,*)(sigii(i),i=1,ntyp)
1098       read (isidep,*)(chip(i),i=1,ntyp)
1099       read (isidep,*)(alp(i),i=1,ntyp)
1100       do i=1,ntyp
1101        read (isidep,*)(epslip(i,j),j=i,ntyp)
1102 C       write(iout,*) "WARNING!!",i,ntyp
1103        write(iout,*) "epslip", i, (epslip(i,j),j=i,ntyp)
1104 C       do j=1,ntyp
1105 C       epslip(i,j)=epslip(i,j)+0.05d0
1106 C       enddo
1107       enddo
1108 C For the GB potential convert sigma'**2 into chi'
1109       if (ipot.eq.4) then
1110         do i=1,ntyp
1111           chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1112         enddo
1113       endif
1114       if (lprint) then
1115         write (iout,'(/a/)') 'Parameters of the BP potential:'
1116         write (iout,'(a/)') 'The epsilon array:'
1117         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1118         write (iout,'(/a)') 'One-body parameters:'
1119         write (iout,'(a,4x,4a)') 'residue','   sigma  ','s||/s_|_^2',
1120      &       '    chip  ','    alph  '
1121         write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1122      &                     chip(i),alp(i),i=1,ntyp)
1123       endif
1124       goto 50
1125 C--------------------- GBV potential -----------------------------------
1126    40 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
1127      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1128      &  (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1129       if (lprint) then
1130         write (iout,'(/a/)') 'Parameters of the GBV potential:'
1131         write (iout,'(a/)') 'The epsilon array:'
1132         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1133         write (iout,'(/a)') 'One-body parameters:'
1134         write (iout,'(a,4x,5a)') 'residue','   sigma  ','    r0    ',
1135      &      's||/s_|_^2','    chip  ','    alph  '
1136         write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1137      &           sigii(i),chip(i),alp(i),i=1,ntyp)
1138       endif
1139    50 continue
1140       close (isidep)
1141 C-----------------------------------------------------------------------
1142 C Calculate the "working" parameters of SC interactions.
1143       do i=2,ntyp
1144         do j=1,i-1
1145           eps(i,j)=eps(j,i)
1146           epslip(i,j)=epslip(j,i)
1147         enddo
1148       enddo
1149       do i=1,ntyp
1150         do j=i,ntyp
1151           sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1152           sigma(j,i)=sigma(i,j)
1153           rs0(i,j)=dwa16*sigma(i,j)
1154           rs0(j,i)=rs0(i,j)
1155         enddo
1156       enddo
1157       if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1158      & 'Working parameters of the SC interactions:',
1159      & '     a    ','     b    ','   augm   ','  sigma ','   r0   ',
1160      & '  chi1   ','   chi2   '
1161       do i=1,ntyp
1162         do j=i,ntyp
1163           epsij=eps(i,j)
1164           epsijlip=epslip(i,j)
1165           if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1166             rrij=sigma(i,j)
1167           else
1168             rrij=rr0(i)+rr0(j)
1169           endif
1170           r0(i,j)=rrij
1171           r0(j,i)=rrij
1172           rrij=rrij**expon
1173           epsij=eps(i,j)
1174           sigeps=dsign(1.0D0,epsij)
1175           epsij=dabs(epsij)
1176           aa_aq(i,j)=epsij*rrij*rrij
1177           bb_aq(i,j)=-sigeps*epsij*rrij
1178           aa_aq(j,i)=aa_aq(i,j)
1179           bb_aq(j,i)=bb_aq(i,j)
1180           sigeps=dsign(1.0D0,epsijlip)
1181           epsijlip=dabs(epsijlip)
1182           aa_lip(i,j)=epsijlip*rrij*rrij
1183           bb_lip(i,j)=-sigeps*epsijlip*rrij
1184           aa_lip(j,i)=aa_lip(i,j)
1185           bb_lip(j,i)=bb_lip(i,j)
1186           if (ipot.gt.2) then
1187             sigt1sq=sigma0(i)**2
1188             sigt2sq=sigma0(j)**2
1189             sigii1=sigii(i)
1190             sigii2=sigii(j)
1191             ratsig1=sigt2sq/sigt1sq
1192             ratsig2=1.0D0/ratsig1
1193             chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1194             if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1195             rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1196           else
1197             rsum_max=sigma(i,j)
1198           endif
1199 c         if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1200             sigmaii(i,j)=rsum_max
1201             sigmaii(j,i)=rsum_max
1202 c         else
1203 c           sigmaii(i,j)=r0(i,j)
1204 c           sigmaii(j,i)=r0(i,j)
1205 c         endif
1206 cd        write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1207           if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1208             r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1209             augm(i,j)=epsij*r_augm**(2*expon)
1210 c           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1211             augm(j,i)=augm(i,j)
1212           else
1213             augm(i,j)=0.0D0
1214             augm(j,i)=0.0D0
1215           endif
1216           if (lprint) then
1217             write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1218      &      restyp(i),restyp(j),aa_aq(i,j),bb_aq(i,j),augm(i,j),
1219      &      sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1220           endif
1221         enddo
1222       enddo
1223 C
1224 C Define the SC-p interaction constants
1225 C
1226 #ifdef OLDSCP
1227       do i=1,20
1228 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1229 C helix formation)
1230 c       aad(i,1)=0.3D0*4.0D0**12
1231 C Following line for constants currently implemented
1232 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1233         aad(i,1)=1.5D0*4.0D0**12
1234 c       aad(i,1)=0.17D0*5.6D0**12
1235         aad(i,2)=aad(i,1)
1236 C "Soft" SC-p repulsion
1237         bad(i,1)=0.0D0
1238 C Following line for constants currently implemented
1239 c       aad(i,1)=0.3D0*4.0D0**6
1240 C "Hard" SC-p repulsion
1241         bad(i,1)=3.0D0*4.0D0**6
1242 c       bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1243         bad(i,2)=bad(i,1)
1244 c       aad(i,1)=0.0D0
1245 c       aad(i,2)=0.0D0
1246 c       bad(i,1)=1228.8D0
1247 c       bad(i,2)=1228.8D0
1248       enddo
1249 #else
1250 C
1251 C 8/9/01 Read the SC-p interaction constants from file
1252 C
1253       do i=1,ntyp
1254         read (iscpp,*) (eps_scp(i,j),rscp(i,j),j=1,2)
1255       enddo
1256       do i=1,ntyp
1257         aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1258         aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1259         bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1260         bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1261       enddo
1262
1263       if (lprint) then
1264         write (iout,*) "Parameters of SC-p interactions:"
1265         do i=1,20
1266           write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1267      &     eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1268         enddo
1269       endif
1270 #endif
1271 C
1272 C Define the constants of the disulfide bridge
1273 C
1274 C      ebr=-12.0D0
1275 c
1276 c Old arbitrary potential - commented out.
1277 c
1278 c      dbr= 4.20D0
1279 c      fbr= 3.30D0
1280 c
1281 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1282 c energy surface of diethyl disulfide.
1283 c A. Liwo and U. Kozlowska, 11/24/03
1284 c
1285 C      D0CM = 3.78d0
1286 C      AKCM = 15.1d0
1287 C      AKTH = 11.0d0
1288 C      AKCT = 12.0d0
1289 C      V1SS =-1.08d0
1290 C      V2SS = 7.61d0
1291 C      V3SS = 13.7d0
1292       write (iout,*) dyn_ss,'dyndyn'
1293       if (dyn_ss) then
1294         ss_depth=ebr/wsc-0.25*eps(1,1)
1295 C        write(iout,*) akcm,whpb,wsc,'KURWA'
1296         Ht=Ht/wsc-0.25*eps(1,1)
1297
1298         akcm=akcm*whpb/wsc
1299         akth=akth*whpb/wsc
1300         akct=akct*whpb/wsc
1301         v1ss=v1ss*whpb/wsc
1302         v2ss=v2ss*whpb/wsc
1303         v3ss=v3ss*whpb/wsc
1304       else
1305         ss_depth=ebr/whpb-0.25*eps(1,1)*wsc/whpb
1306       endif
1307
1308 C      if (lprint) then
1309       write (iout,'(/a)') "Disulfide bridge parameters:"
1310       write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1311       write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1312       write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1313       write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1314      & ' v3ss:',v3ss
1315 C      endif
1316       return
1317       end