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