source/maxlik/src-Fmatch/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 'COMMON.IOUNITS'
  10       include 'COMMON.CHAIN'
  11       include 'COMMON.INTERACT'
  12       include 'COMMON.GEO'
  13       include 'COMMON.LOCAL'
  14       include 'COMMON.TORSION'
  15       include 'COMMON.FFIELD'
  16       include 'COMMON.NAMES'
  17       include 'COMMON.SBRIDGE'
  18       include 'COMMON.WEIGHTS'
  19       include 'COMMON.SCCOR'
  20       include 'COMMON.SCROT'
  21       include 'COMMON.SHIELD'
  22       include 'COMMON.CONTROL'
  23       include 'COMMON.DERIV'
  24       include 'COMMON.OPTIM'
  25       include 'COMMON.ENERGIES'
  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       double precision ip,mp
  40       character*6 res1
  41       character*3 lancuch,ucase
  42 C
  43 C Body
  44 C
  45       call getenv("PRINT_PARM",lancuch)
  46       lprint = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
  47       write (iout,*) "lprint ",lprint
  48 C Set LPRINT=.TRUE. for debugging
  49       dwa16=2.0d0**(1.0d0/6.0d0)
  50       itypro=20
  51 C Assign virtual-bond length
  52       vbl=3.8D0
  53       vblinv=1.0D0/vbl
  54       vblinv2=vblinv*vblinv
  55       call card_concat(controlcard,.true.)
  56       wname(4)="WCORRH"
  57       do i=1,n_ene
  58         key = wname(i)(:ilen(wname(i)))
  59         call reada(controlcard,key(:ilen(key)),ww(i),0.0d0)
  60       enddo
  61 c      print *,me," checkpoint 7"
  62
  63       write (iout,*) "iparm",iparm," myparm",myparm
  64       call flush(iout)
  65 c If reading not own parameters, skip assignment
  66       call reada(controlcard,"D0CM",d0cm,3.78d0)
  67       call reada(controlcard,"AKCM",akcm,15.1d0)
  68       call reada(controlcard,"AKTH",akth,11.0d0)
  69       call reada(controlcard,"AKCT",akct,12.0d0)
  70       call reada(controlcard,"V1SS",v1ss,-1.08d0)
  71       call reada(controlcard,"V2SS",v2ss,7.61d0)
  72       call reada(controlcard,"V3SS",v3ss,13.7d0)
  73       call reada(controlcard,"EBR",ebr,-5.50D0)
  74       call reada(controlcard,"DTRISS",dtriss,1.0D0)
  75       call reada(controlcard,"ATRISS",atriss,0.3D0)
  76       call reada(controlcard,"BTRISS",btriss,0.02D0)
  77       call reada(controlcard,"CTRISS",ctriss,1.0D0)
  78       dyn_ss=(index(controlcard,'DYN_SS').gt.0)
  79       write(iout,*) "ATRISS",atriss
  80       write(iout,*) "BTRISS",btriss
  81       write(iout,*) "CTRISS",ctriss
  82       write(iout,*) "DTRISS",dtriss
  83
  84 C      do i=1,maxres
  85 C        dyn_ss_mask(i)=.false.
  86 C      enddo
  87 C      ebr=-12.0D0
  88 c
  89 c Old arbitrary potential - commented out.
  90 c
  91 c      dbr= 4.20D0
  92 c      fbr= 3.30D0
  93 c
  94 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
  95 c energy surface of diethyl disulfide.
  96 c A. Liwo and U. Kozlowska, 11/24/03
  97 c
  98       D0CM = 3.78d0
  99       AKCM = 15.1d0
 100       AKTH = 11.0d0
 101       AKCT = 12.0d0
 102       V1SS =-1.08d0
 103       V2SS = 7.61d0
 104       V3SS = 13.7d0
 105
 106       do i=1,maxres-1
 107         do j=i+1,maxres
 108           dyn_ssbond_ij(i,j)=1.0d300
 109         enddo
 110       enddo
 111       call reada(controlcard,"HT",Ht,0.0D0)
 112 C      if (dyn_ss) then
 113 C        ss_depth=ebr/wsc-0.25*eps(1,1)
 114 C        write(iout,*) HT,wsc,eps(1,1),'KURWA'
 115 C        Ht=Ht/wsc-0.25*eps(1,1)
 116
 117 C        akcm=akcm*whpb/wsc
 118 C        akth=akth*whpb/wsc
 119 C        akct=akct*whpb/wsc
 120 C        v1ss=v1ss*whpb/wsc
 121 C        v2ss=v2ss*whpb/wsc
 122 C        v3ss=v3ss*whpb/wsc
 123 C      else
 124 C        ss_depth=ebr/whpb-0.25*eps(1,1)*wsc/whpb
 125 C      endif
 126
 127       if (iparm.eq.myparm) then
 128
 129 c
 130 c Setup weights for UNRES
 131 c
 132       wsc=ww(1)
 133       wscp=ww(2)
 134       welec=ww(3)
 135       wcorr=ww(4)
 136       wcorr5=ww(5)
 137       wcorr6=ww(6)
 138       wel_loc=ww(7)
 139       wturn3=ww(8)
 140       wturn4=ww(9)
 141       wturn6=ww(10)
 142       wang=ww(11)
 143       wscloc=ww(12)
 144       wtor=ww(13)
 145       wtor_d=ww(14)
 146       wvdwpp=ww(16)
 147       wbond=ww(18)
 148       wsccor=ww(19)
 149       whpb=ww(15)
 150       wstrain=ww(15)
 151       wliptran=ww(22)
 152       wshield=ww(25)
 153       endif
 154
 155       call card_concat(controlcard,.false.)
 156
 157 c Return if not own parameters
 158
 159       if (iparm.ne.myparm) return
 160
 161       call reads(controlcard,"BONDPAR",bondname_t,bondname)
 162       open (ibond,file=bondname_t,status='old')
 163       rewind(ibond)
 164       call reads(controlcard,"THETPAR",thetname_t,thetname)
 165       open (ithep,file=thetname_t,status='old')
 166       rewind(ithep)
 167       call reads(controlcard,"ROTPAR",rotname_t,rotname)
 168       open (irotam,file=rotname_t,status='old')
 169       rewind(irotam)
 170       call reads(controlcard,"TORPAR",torname_t,torname)
 171       open (itorp,file=torname_t,status='old')
 172       rewind(itorp)
 173       call reads(controlcard,"TORDPAR",tordname_t,tordname)
 174       write (iout,*) "tor_mode",tor_mode
 175       call flush(iout)
 176       if (tor_mode.eq.0)
 177      & open (itordp,file=tordname_t,status='old')
 178       rewind(itordp)
 179       call reads(controlcard,"SCCORPAR",sccorname_t,sccorname)
 180       open (isccor,file=sccorname_t,status='old')
 181       rewind(isccor)
 182       call reads(controlcard,"FOURIER",fouriername_t,fouriername)
 183       open (ifourier,file=fouriername_t,status='old')
 184       rewind(ifourier)
 185       call reads(controlcard,"ELEPAR",elename_t,elename)
 186       open (ielep,file=elename_t,status='old')
 187       rewind(ielep)
 188       call reads(controlcard,"SIDEPAR",sidename_t,sidename)
 189       open (isidep,file=sidename_t,status='old')
 190       rewind(isidep)
 191       call reads(controlcard,"SCPPAR",scpname_t,scpname)
 192       open (iscpp,file=scpname_t,status='old')
 193       rewind(iscpp)
 194       write (iout,*) "Parameter set:",iparm
 195       write (iout,*) "Energy-term weights:"
 196       do i=1,n_ene
 197         write (iout,'(a16,f10.5)') wname(i),ww(i)
 198       enddo
 199       write (iout,*) "Sidechain potential file        : ",
 200      &  sidename_t(:ilen(sidename_t))
 201 #ifndef OLDSCP
 202       write (iout,*) "SCp potential file              : ",
 203      &  scpname_t(:ilen(scpname_t))
 204 #endif
 205       write (iout,*) "Electrostatic potential file    : ",
 206      &  elename_t(:ilen(elename_t))
 207       write (iout,*) "Cumulant coefficient file       : ",
 208      &  fouriername_t(:ilen(fouriername_t))
 209       write (iout,*) "Torsional parameter file        : ",
 210      &  torname_t(:ilen(torname_t))
 211       write (iout,*) "Double torsional parameter file : ",
 212      &  tordname_t(:ilen(tordname_t))
 213       write (iout,*) "Backbone-rotamer parameter file : ",
 214      &  sccorname(:ilen(sccorname))
 215       write (iout,*) "Bond & inertia constant file    : ",
 216      &  bondname_t(:ilen(bondname_t))
 217       write (iout,*) "Bending parameter file          : ",
 218      &  thetname_t(:ilen(thetname_t))
 219       write (iout,*) "Rotamer parameter file          : ",
 220      &  rotname_t(:ilen(rotname_t))
 221
 222 c
 223 c Read the virtual-bond parameters, masses, and moments of inertia
 224 c and Stokes' radii of the peptide group and side chains
 225 c
 226 #ifdef CRYST_BOND
 227       read (ibond,*,end=121,err=121) vbldp0,vbldpdum,akp
 228       do i=1,ntyp
 229         nbondterm(i)=1
 230         read (ibond,*,end=121,err=121) vbldsc0(1,i),aksc(1,i)
 231         dsc(i) = vbldsc0(1,i)
 232         if (i.eq.10) then
 233           dsc_inv(i)=0.0D0
 234         else
 235           dsc_inv(i)=1.0D0/dsc(i)
 236         endif
 237       enddo
 238 #else
 239       read (ibond,*,end=121,err=121) ijunk,vbldp0,vbldpdum,akp,rjunk
 240       do i=1,ntyp
 241         read (ibond,*,end=121,err=121) nbondterm(i),(vbldsc0(j,i),
 242      &   aksc(j,i),abond0(j,i),j=1,nbondterm(i))
 243         dsc(i) = vbldsc0(1,i)
 244         if (i.eq.10) then
 245           dsc_inv(i)=0.0D0
 246         else
 247           dsc_inv(i)=1.0D0/dsc(i)
 248         endif
 249       enddo
 250 #endif
 251       if (lprint) then
 252         write(iout,'(/a/)')"Force constants virtual bonds:"
 253         write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K',
 254      &   'inertia','Pstok'
 255         write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0
 256         do i=1,ntyp
 257           write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),
 258      &      vbldsc0(1,i),aksc(1,i),abond0(1,i)
 259           do j=2,nbondterm(i)
 260             write (iout,'(13x,3f10.5)')
 261      &        vbldsc0(j,i),aksc(j,i),abond0(j,i)
 262           enddo
 263         enddo
 264       endif
 265        read(iliptranpar,*,end=1161,err=1161) pepliptran
 266        do i=1,ntyp
 267        read(iliptranpar,*,end=1161,err=1161) liptranene(i)
 268        enddo
 269        rewind(iliptranpar)
 270 #ifdef CRYST_THETA
 271 C
 272 C Read the parameters of the probability distribution/energy expression
 273 C of the virtual-bond valence angles theta
 274 C
 275       do i=1,ntyp
 276         read (ithep,*,end=111,err=111) a0thet(i),(athet(j,i,1,1),j=1,2),
 277      &    (bthet(j,i,1,1),j=1,2)
 278         read (ithep,*,end=111,err=111) (polthet(j,i),j=0,3)
 279         read (ithep,*,end=111,err=111) (gthet(j,i),j=1,3)
 280         read (ithep,*,end=111,err=111) theta0(i),sig0(i),sigc0(i)
 281         sigc0(i)=sigc0(i)**2
 282       enddo
 283       do i=1,ntyp
 284       athet(1,i,1,-1)=athet(1,i,1,1)
 285       athet(2,i,1,-1)=athet(2,i,1,1)
 286       bthet(1,i,1,-1)=-bthet(1,i,1,1)
 287       bthet(2,i,1,-1)=-bthet(2,i,1,1)
 288       athet(1,i,-1,1)=-athet(1,i,1,1)
 289       athet(2,i,-1,1)=-athet(2,i,1,1)
 290       bthet(1,i,-1,1)=bthet(1,i,1,1)
 291       bthet(2,i,-1,1)=bthet(2,i,1,1)
 292       enddo
 293       do i=-ntyp,-1
 294       a0thet(i)=a0thet(-i)
 295       athet(1,i,-1,-1)=athet(1,-i,1,1)
 296       athet(2,i,-1,-1)=-athet(2,-i,1,1)
 297       bthet(1,i,-1,-1)=bthet(1,-i,1,1)
 298       bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
 299       athet(1,i,-1,1)=athet(1,-i,1,1)
 300       athet(2,i,-1,1)=-athet(2,-i,1,1)
 301       bthet(1,i,-1,1)=-bthet(1,-i,1,1)
 302       bthet(2,i,-1,1)=bthet(2,-i,1,1)
 303       athet(1,i,1,-1)=-athet(1,-i,1,1)
 304       athet(2,i,1,-1)=athet(2,-i,1,1)
 305       bthet(1,i,1,-1)=bthet(1,-i,1,1)
 306       bthet(2,i,1,-1)=-bthet(2,-i,1,1)
 307       theta0(i)=theta0(-i)
 308       sig0(i)=sig0(-i)
 309       sigc0(i)=sigc0(-i)
 310        do j=0,3
 311         polthet(j,i)=polthet(j,-i)
 312        enddo
 313        do j=1,3
 314          gthet(j,i)=gthet(j,-i)
 315        enddo
 316       enddo
 317       close (ithep)
 318       if (lprint) then
 319 c       write (iout,'(a)')
 320 c    &    'Parameters of the virtual-bond valence angles:'
 321 c       write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
 322 c    & '    ATHETA0   ','         A1   ','        A2    ',
 323 c    & '        B1    ','         B2   '
 324 c       do i=1,ntyp
 325 c         write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
 326 c    &        a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2)
 327 c       enddo
 328 c       write (iout,'(/a/9x,5a/79(1h-))')
 329 c    & 'Parameters of the expression for sigma(theta_c):',
 330 c    & '     ALPH0    ','      ALPH1   ','     ALPH2    ',
 331 c    & '     ALPH3    ','    SIGMA0C   '
 332 c       do i=1,ntyp
 333 c         write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
 334 c    &      (polthet(j,i),j=0,3),sigc0(i)
 335 c       enddo
 336 c       write (iout,'(/a/9x,5a/79(1h-))')
 337 c    & 'Parameters of the second gaussian:',
 338 c    & '    THETA0    ','     SIGMA0   ','        G1    ',
 339 c    & '        G2    ','         G3   '
 340 c       do i=1,ntyp
 341 c         write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
 342 c    &       sig0(i),(gthet(j,i),j=1,3)
 343 c       enddo
 344         write (iout,'(a)')
 345      &    'Parameters of the virtual-bond valence angles:'
 346         write (iout,'(/a/9x,5a/79(1h-))')
 347      & 'Coefficients of expansion',
 348      & '     theta0   ','    a1*10^2   ','   a2*10^2    ',
 349      & '   b1*10^1    ','    b2*10^1   '
 350         do i=1,ntyp
 351           write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),
 352      &        a0thet(i),(100*athet(j,i,1,1),j=1,2),
 353      &        (10*bthet(j,i,1,1),j=1,2)
 354         enddo
 355         write (iout,'(/a/9x,5a/79(1h-))')
 356      & 'Parameters of the expression for sigma(theta_c):',
 357      & ' alpha0       ','  alph1       ',' alph2        ',
 358      & ' alhp3        ','   sigma0c    '
 359         do i=1,ntyp
 360           write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),
 361      &      (polthet(j,i),j=0,3),sigc0(i)
 362         enddo
 363         write (iout,'(/a/9x,5a/79(1h-))')
 364      & 'Parameters of the second gaussian:',
 365      & '    theta0    ','  sigma0*10^2 ','      G1*10^-1',
 366      & '        G2    ','   G3*10^1    '
 367         do i=1,ntyp
 368           write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),
 369      &       100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
 370         enddo
 371       endif
 372 #else
 373       IF (tor_mode.eq.0) THEN
 374 C
 375 C Read the parameters of Utheta determined from ab initio surfaces
 376 C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
 377 C
 378       read (ithep,*,end=111,err=111) nthetyp,ntheterm,ntheterm2,
 379      &  ntheterm3,nsingle,ndouble
 380       nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
 381       read (ithep,*,end=111,err=111) (ithetyp(i),i=1,ntyp1)
 382       do i=-ntyp1,-1
 383         ithetyp(i)=-ithetyp(-i)
 384       enddo
 385       do iblock=1,2
 386       do i=-maxthetyp,maxthetyp
 387         do j=-maxthetyp,maxthetyp
 388           do k=-maxthetyp,maxthetyp
 389             aa0thet(i,j,k,iblock)=0.0d0
 390             do l=1,ntheterm
 391               aathet(l,i,j,k,iblock)=0.0d0
 392             enddo
 393             do l=1,ntheterm2
 394               do m=1,nsingle
 395                 bbthet(m,l,i,j,k,iblock)=0.0d0
 396                 ccthet(m,l,i,j,k,iblock)=0.0d0
 397                 ddthet(m,l,i,j,k,iblock)=0.0d0
 398                 eethet(m,l,i,j,k,iblock)=0.0d0
 399               enddo
 400             enddo
 401             do l=1,ntheterm3
 402               do m=1,ndouble
 403                 do mm=1,ndouble
 404                  ffthet(mm,m,l,i,j,k,iblock)=0.0d0
 405                  ggthet(mm,m,l,i,j,k,iblock)=0.0d0
 406                 enddo
 407               enddo
 408             enddo
 409           enddo
 410         enddo
 411       enddo
 412       enddo
 413 C      write (iout,*) "KURWA1"
 414       do iblock=1,2
 415       do i=0,nthetyp
 416         do j=-nthetyp,nthetyp
 417           do k=-nthetyp,nthetyp
 418             read (ithep,'(6a)',end=111,err=111) res1
 419             write(iout,*) res1,i,j,k
 420             read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock)
 421             read (ithep,*,end=111,err=111)(aathet(l,i,j,k,iblock),
 422      &        l=1,ntheterm)
 423             read (ithep,*,end=111,err=111)
 424      &       ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
 425      &        (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
 426      &        (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
 427      &        (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle)
 428      &        ,ll=1,ntheterm2)
 429             read (ithep,*,end=111,err=111)
 430      &      (((ffthet(llll,lll,ll,i,j,k,iblock),
 431      &      ffthet(lll,llll,ll,i,j,k,iblock),
 432      &         ggthet(llll,lll,ll,i,j,k,iblock)
 433      &        ,ggthet(lll,llll,ll,i,j,k,iblock),
 434      &         llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
 435           enddo
 436         enddo
 437       enddo
 438 C       write(iout,*) "KURWA1.1"
 439 C
 440 C For dummy ends assign glycine-type coefficients of theta-only terms; the
 441 C coefficients of theta-and-gamma-dependent terms are zero.
 442 C
 443       do i=1,nthetyp
 444         do j=1,nthetyp
 445           do l=1,ntheterm
 446             aathet(l,i,j,nthetyp+1,iblock)=0.0d0
 447             aathet(l,nthetyp+1,i,j,iblock)=0.0d0
 448           enddo
 449           aa0thet(i,j,nthetyp+1,iblock)=0.0d0
 450           aa0thet(nthetyp+1,i,j,iblock)=0.0d0
 451         enddo
 452         do l=1,ntheterm
 453           aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
 454         enddo
 455         aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
 456       enddo
 457       enddo
 458 C       write(iout,*) "KURWA1.5"
 459 C Substitution for D aminoacids from symmetry.
 460       do iblock=1,2
 461       do i=-nthetyp,0
 462         do j=-nthetyp,nthetyp
 463           do k=-nthetyp,nthetyp
 464            aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
 465            do l=1,ntheterm
 466            aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
 467            enddo
 468            do ll=1,ntheterm2
 469             do lll=1,nsingle
 470             bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
 471             ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
 472             ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
 473             eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
 474             enddo
 475           enddo
 476           do ll=1,ntheterm3
 477            do lll=2,ndouble
 478             do llll=1,lll-1
 479             ffthet(llll,lll,ll,i,j,k,iblock)=
 480      &      ffthet(llll,lll,ll,-i,-j,-k,iblock)
 481             ffthet(lll,llll,ll,i,j,k,iblock)=
 482      &      ffthet(lll,llll,ll,-i,-j,-k,iblock)
 483             ggthet(llll,lll,ll,i,j,k,iblock)=
 484      &      -ggthet(llll,lll,ll,-i,-j,-k,iblock)
 485             ggthet(lll,llll,ll,i,j,k,iblock)=
 486      &      -ggthet(lll,llll,ll,-i,-j,-k,iblock)
 487             enddo !ll
 488            enddo  !lll
 489           enddo   !llll
 490          enddo    !k
 491         enddo     !j
 492        enddo      !i
 493       enddo       !iblock
 494
 495 C
 496 C Control printout of the coefficients of virtual-bond-angle potentials
 497 C
 498       if (lprint) then
 499         write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
 500         do iblock=1,2
 501         do i=1,nthetyp+1
 502           do j=1,nthetyp+1
 503             do k=1,nthetyp+1
 504               write (iout,'(//4a)')
 505      &         'Type ',onelett(i),onelett(j),onelett(k)
 506               write (iout,'(//a,10x,a)') " l","a[l]"
 507               write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
 508               write (iout,'(i2,1pe15.5)')
 509      &           (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
 510             do l=1,ntheterm2
 511               write (iout,'(//2h m,4(9x,a,3h[m,i1,1h]))')
 512      &          "b",l,"c",l,"d",l,"e",l
 513               do m=1,nsingle
 514                 write (iout,'(i2,4(1pe15.5))') m,
 515      &          bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),
 516      &          ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
 517               enddo
 518             enddo
 519             do l=1,ntheterm3
 520               write (iout,'(//3hm,n,4(6x,a,5h[m,n,i1,1h]))')
 521      &          "f+",l,"f-",l,"g+",l,"g-",l
 522               do m=2,ndouble
 523                 do n=1,m-1
 524                   write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
 525      &              ffthet(n,m,l,i,j,k,iblock),
 526      &              ffthet(m,n,l,i,j,k,iblock),
 527      &              ggthet(n,m,l,i,j,k,iblock),
 528      &              ggthet(m,n,l,i,j,k,iblock)
 529                 enddo
 530               enddo
 531             enddo
 532           enddo
 533         enddo
 534       enddo
 535       enddo
 536       call flush(iout)
 537       endif
 538
 539       ELSE
 540
 541 C here will be the apropriate recalibrating for D-aminoacid
 542       read (ithep,*,end=111,err=111) nthetyp
 543       do i=-nthetyp+1,nthetyp-1
 544         read (ithep,*,end=111,err=111) nbend_kcc_Tb(i)
 545         do j=0,nbend_kcc_Tb(i)
 546           read (ithep,*,end=111,err=111) ijunk,v1bend_chyb(j,i)
 547         enddo
 548       enddo
 549       if (lprint) then
 550         write (iout,'(a)')
 551      &    "Parameters of the valence-only potentials"
 552         do i=-nthetyp+1,nthetyp-1
 553         write (iout,'(2a)') "Type ",toronelet(i)
 554         do k=0,nbend_kcc_Tb(i)
 555           write(iout,'(i5,f15.5)') k,v1bend_chyb(k,i)
 556         enddo
 557         enddo
 558       endif
 559
 560       ENDIF ! TOR_MODE
 561
 562       close(ithep)
 563 #endif
 564 C      write(iout,*) 'KURWA2'
 565 #ifdef CRYST_SC
 566 C
 567 C Read the parameters of the probability distribution/energy expression
 568 C of the side chains.
 569 C
 570       do i=1,ntyp
 571 cc      write (iout,*) "tu dochodze",i
 572         read (irotam,'(3x,i3,f8.3)') nlob(i),dsc(i)
 573         if (i.eq.10) then
 574           dsc_inv(i)=0.0D0
 575         else
 576           dsc_inv(i)=1.0D0/dsc(i)
 577         endif
 578         if (i.ne.10) then
 579         do j=1,nlob(i)
 580           do k=1,3
 581             do l=1,3
 582               blower(l,k,j)=0.0D0
 583             enddo
 584           enddo
 585         enddo
 586         bsc(1,i)=0.0D0
 587         read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),
 588      &    ((blower(k,l,1),l=1,k),k=1,3)
 589         censc(1,1,-i)=censc(1,1,i)
 590         censc(2,1,-i)=censc(2,1,i)
 591         censc(3,1,-i)=-censc(3,1,i)
 592         do j=2,nlob(i)
 593           read (irotam,*,end=112,err=112) bsc(j,i)
 594           read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),
 595      &                                 ((blower(k,l,j),l=1,k),k=1,3)
 596         censc(1,j,-i)=censc(1,j,i)
 597         censc(2,j,-i)=censc(2,j,i)
 598         censc(3,j,-i)=-censc(3,j,i)
 599 C BSC is amplitude of Gaussian
 600         enddo
 601         do j=1,nlob(i)
 602           do k=1,3
 603             do l=1,k
 604               akl=0.0D0
 605               do m=1,3
 606                 akl=akl+blower(k,m,j)*blower(l,m,j)
 607               enddo
 608               gaussc(k,l,j,i)=akl
 609               gaussc(l,k,j,i)=akl
 610              if (((k.eq.3).and.(l.ne.3))
 611      &        .or.((l.eq.3).and.(k.ne.3))) then
 612                 gaussc(k,l,j,-i)=-akl
 613                 gaussc(l,k,j,-i)=-akl
 614               else
 615                 gaussc(k,l,j,-i)=akl
 616                 gaussc(l,k,j,-i)=akl
 617               endif
 618             enddo
 619           enddo
 620         enddo
 621         endif
 622       enddo
 623       close (irotam)
 624       if (lprint) then
 625         write (iout,'(/a)') 'Parameters of side-chain local geometry'
 626         do i=1,ntyp
 627           nlobi=nlob(i)
 628           if (nlobi.gt.0) then
 629           write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
 630      &     ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
 631 c          write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
 632 c          write (iout,'(a,f10.4,4(16x,f10.4))')
 633 c     &                             'Center  ',(bsc(j,i),j=1,nlobi)
 634 c          write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),j=1,nlobi)
 635            write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
 636      &                             'log h',(bsc(j,i),j=1,nlobi)
 637            write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
 638      &    'x',((censc(k,j,i),k=1,3),j=1,nlobi)
 639 c          write (iout,'(a)')
 640 c         do j=1,nlobi
 641 c           ind=0
 642 c           do k=1,3
 643 c             do l=1,k
 644 c              ind=ind+1
 645 c              blower(k,l,j)=gaussc(ind,j,i)
 646 c             enddo
 647 c           enddo
 648 c         enddo
 649           do k=1,3
 650             write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
 651      &                 ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
 652           enddo
 653           endif
 654         enddo
 655       endif
 656 #else
 657 C
 658 C Read scrot parameters for potentials determined from all-atom AM1 calculations
 659 C added by Urszula Kozlowska 07/11/2007
 660 C
 661       do i=1,ntyp
 662         read (irotam,*,end=112,err=112)
 663        if (i.eq.10) then
 664          read (irotam,*,end=112,err=112)
 665        else
 666          do j=1,65
 667            read(irotam,*,end=112,err=112) sc_parmin(j,i)
 668          enddo
 669        endif
 670       enddo
 671 #endif
 672       close(irotam)
 673 C
 674 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
 675 C         interaction energy of the Gly, Ala, and Pro prototypes.
 676 C
 677       read (ifourier,*,end=115,err=115) nloctyp
 678       SPLIT_FOURIERTOR = nloctyp.lt.0
 679       nloctyp = iabs(nloctyp)
 680       write (iout,*) "nloctyp",nloctyp
 681 #ifdef NEWCORR
 682       read (ifourier,*,end=115,err=115) (itype2loc(i),i=1,ntyp)
 683       read (ifourier,*,end=115,err=115) (iloctyp(i),i=0,nloctyp-1)
 684       itype2loc(ntyp1)=nloctyp
 685       iloctyp(nloctyp)=ntyp1
 686       do i=1,ntyp1
 687         itype2loc(-i)=-itype2loc(i)
 688       enddo
 689 #else
 690       iloctyp(0)=10
 691       iloctyp(1)=9
 692       iloctyp(2)=20
 693       iloctyp(3)=ntyp1
 694 #endif
 695       do i=1,nloctyp
 696         iloctyp(-i)=-iloctyp(i)
 697       enddo
 698 c      write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1)
 699 c      write (iout,*) "nloctyp",nloctyp,
 700 c     &  " iloctyp",(iloctyp(i),i=0,nloctyp)
 701 #ifdef NEWCORR
 702       do i=0,nloctyp-1
 703 c             write (iout,*) "NEWCORR",i
 704         read (ifourier,*,end=115,err=115)
 705         do ii=1,3
 706           do j=1,2
 707             read (ifourier,*,end=115,err=115) bnew1(ii,j,i)
 708           enddo
 709         enddo
 710 c             write (iout,*) "NEWCORR BNEW1"
 711 c             write (iout,*) ((bnew1(ii,j,i),ii=1,3),j=1,2)
 712         do ii=1,3
 713           do j=1,2
 714             read (ifourier,*,end=115,err=115) bnew2(ii,j,i)
 715           enddo
 716         enddo
 717 c             write (iout,*) "NEWCORR BNEW2"
 718 c             write (iout,*) ((bnew2(ii,j,i),ii=1,3),j=1,2)
 719         do kk=1,3
 720           read (ifourier,*,end=115,err=115) ccnew(kk,1,i)
 721           read (ifourier,*,end=115,err=115) ccnew(kk,2,i)
 722         enddo
 723 c             write (iout,*) "NEWCORR CCNEW"
 724 c             write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
 725         do kk=1,3
 726           read (ifourier,*,end=115,err=115) ddnew(kk,1,i)
 727           read (ifourier,*,end=115,err=115) ddnew(kk,2,i)
 728         enddo
 729 c             write (iout,*) "NEWCORR DDNEW"
 730 c             write (iout,*) ((ddnew(ii,j,i),ii=1,3),j=1,2)
 731         do ii=1,2
 732           do jj=1,2
 733             do kk=1,2
 734               read (ifourier,*,end=115,err=115) eenew(ii,jj,kk,i)
 735             enddo
 736           enddo
 737         enddo
 738 c             write (iout,*) "NEWCORR EENEW1"
 739 c             write(iout,*)(((eenew(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
 740         do ii=1,3
 741           read (ifourier,*,end=115,err=115) e0new(ii,i)
 742         enddo
 743 c             write (iout,*) (e0new(ii,i),ii=1,3)
 744       enddo
 745 c             write (iout,*) "NEWCORR EENEW"
 746       do i=0,nloctyp-1
 747         do ii=1,3
 748           ccnew(ii,1,i)=ccnew(ii,1,i)/2
 749           ccnew(ii,2,i)=ccnew(ii,2,i)/2
 750           ddnew(ii,1,i)=ddnew(ii,1,i)/2
 751           ddnew(ii,2,i)=ddnew(ii,2,i)/2
 752         enddo
 753       enddo
 754       do i=1,nloctyp-1
 755         do ii=1,3
 756           bnew1(ii,1,-i)= bnew1(ii,1,i)
 757           bnew1(ii,2,-i)=-bnew1(ii,2,i)
 758           bnew2(ii,1,-i)= bnew2(ii,1,i)
 759           bnew2(ii,2,-i)=-bnew2(ii,2,i)
 760         enddo
 761         do ii=1,3
 762 c          ccnew(ii,1,i)=ccnew(ii,1,i)/2
 763 c          ccnew(ii,2,i)=ccnew(ii,2,i)/2
 764 c          ddnew(ii,1,i)=ddnew(ii,1,i)/2
 765 c          ddnew(ii,2,i)=ddnew(ii,2,i)/2
 766           ccnew(ii,1,-i)=ccnew(ii,1,i)
 767           ccnew(ii,2,-i)=-ccnew(ii,2,i)
 768           ddnew(ii,1,-i)=ddnew(ii,1,i)
 769           ddnew(ii,2,-i)=-ddnew(ii,2,i)
 770         enddo
 771         e0new(1,-i)= e0new(1,i)
 772         e0new(2,-i)=-e0new(2,i)
 773         e0new(3,-i)=-e0new(3,i)
 774         do kk=1,2
 775           eenew(kk,1,1,-i)= eenew(kk,1,1,i)
 776           eenew(kk,1,2,-i)=-eenew(kk,1,2,i)
 777           eenew(kk,2,1,-i)=-eenew(kk,2,1,i)
 778           eenew(kk,2,2,-i)= eenew(kk,2,2,i)
 779         enddo
 780       enddo
 781       if (lprint) then
 782         write (iout,'(a)') "Coefficients of the multibody terms"
 783         do i=-nloctyp+1,nloctyp-1
 784           write (iout,*) "Type: ",onelet(iloctyp(i))
 785           write (iout,*) "Coefficients of the expansion of B1"
 786           do j=1,2
 787             write (iout,'(3hB1(,i1,1h),3f10.5)') j,(bnew1(k,j,i),k=1,3)
 788           enddo
 789           write (iout,*) "Coefficients of the expansion of B2"
 790           do j=1,2
 791             write (iout,'(3hB2(,i1,1h),3f10.5)') j,(bnew2(k,j,i),k=1,3)
 792           enddo
 793           write (iout,*) "Coefficients of the expansion of C"
 794           write (iout,'(3hC11,3f10.5)') (ccnew(j,1,i),j=1,3)
 795           write (iout,'(3hC12,3f10.5)') (ccnew(j,2,i),j=1,3)
 796           write (iout,*) "Coefficients of the expansion of D"
 797           write (iout,'(3hD11,3f10.5)') (ddnew(j,1,i),j=1,3)
 798           write (iout,'(3hD12,3f10.5)') (ddnew(j,2,i),j=1,3)
 799           write (iout,*) "Coefficients of the expansion of E"
 800           write (iout,'(2hE0,3f10.5)') (e0new(j,i),j=1,3)
 801           do j=1,2
 802             do k=1,2
 803               write (iout,'(1hE,2i1,2f10.5)') j,k,(eenew(l,j,k,i),l=1,2)
 804             enddo
 805           enddo
 806         enddo
 807       endif
 808       IF (SPLIT_FOURIERTOR) THEN
 809       do i=0,nloctyp-1
 810 c             write (iout,*) "NEWCORR TOR",i
 811         read (ifourier,*,end=115,err=115)
 812         do ii=1,3
 813           do j=1,2
 814             read (ifourier,*,end=115,err=115) bnew1tor(ii,j,i)
 815           enddo
 816         enddo
 817 c             write (iout,*) "NEWCORR BNEW1 TOR"
 818 c             write (iout,*) ((bnew1tor(ii,j,i),ii=1,3),j=1,2)
 819         do ii=1,3
 820           do j=1,2
 821             read (ifourier,*,end=115,err=115) bnew2tor(ii,j,i)
 822           enddo
 823         enddo
 824 c             write (iout,*) "NEWCORR BNEW2 TOR"
 825 c             write (iout,*) ((bnew2tor(ii,j,i),ii=1,3),j=1,2)
 826         do kk=1,3
 827           read (ifourier,*,end=115,err=115) ccnewtor(kk,1,i)
 828           read (ifourier,*,end=115,err=115) ccnewtor(kk,2,i)
 829         enddo
 830 c             write (iout,*) "NEWCORR CCNEW TOR"
 831 c             write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
 832         do kk=1,3
 833           read (ifourier,*,end=115,err=115) ddnewtor(kk,1,i)
 834           read (ifourier,*,end=115,err=115) ddnewtor(kk,2,i)
 835         enddo
 836 c             write (iout,*) "NEWCORR DDNEW TOR"
 837 c             write (iout,*) ((ddnewtor(ii,j,i),ii=1,3),j=1,2)
 838         do ii=1,2
 839           do jj=1,2
 840             do kk=1,2
 841               read (ifourier,*,end=115,err=115) eenewtor(ii,jj,kk,i)
 842             enddo
 843           enddo
 844         enddo
 845 c         write (iout,*) "NEWCORR EENEW1 TOR"
 846 c         write(iout,*)(((eenewtor(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
 847         do ii=1,3
 848           read (ifourier,*,end=115,err=115) e0newtor(ii,i)
 849         enddo
 850 c             write (iout,*) (e0newtor(ii,i),ii=1,3)
 851       enddo
 852 c             write (iout,*) "NEWCORR EENEW TOR"
 853       do i=0,nloctyp-1
 854         do ii=1,3
 855           ccnewtor(ii,1,i)=ccnewtor(ii,1,i)/2
 856           ccnewtor(ii,2,i)=ccnewtor(ii,2,i)/2
 857           ddnewtor(ii,1,i)=ddnewtor(ii,1,i)/2
 858           ddnewtor(ii,2,i)=ddnewtor(ii,2,i)/2
 859         enddo
 860       enddo
 861       do i=1,nloctyp-1
 862         do ii=1,3
 863           bnew1tor(ii,1,-i)= bnew1tor(ii,1,i)
 864           bnew1tor(ii,2,-i)=-bnew1tor(ii,2,i)
 865           bnew2tor(ii,1,-i)= bnew2tor(ii,1,i)
 866           bnew2tor(ii,2,-i)=-bnew2tor(ii,2,i)
 867         enddo
 868         do ii=1,3
 869           ccnewtor(ii,1,-i)=ccnewtor(ii,1,i)
 870           ccnewtor(ii,2,-i)=-ccnewtor(ii,2,i)
 871           ddnewtor(ii,1,-i)=ddnewtor(ii,1,i)
 872           ddnewtor(ii,2,-i)=-ddnewtor(ii,2,i)
 873         enddo
 874         e0newtor(1,-i)= e0newtor(1,i)
 875         e0newtor(2,-i)=-e0newtor(2,i)
 876         e0newtor(3,-i)=-e0newtor(3,i)
 877         do kk=1,2
 878           eenewtor(kk,1,1,-i)= eenewtor(kk,1,1,i)
 879           eenewtor(kk,1,2,-i)=-eenewtor(kk,1,2,i)
 880           eenewtor(kk,2,1,-i)=-eenewtor(kk,2,1,i)
 881           eenewtor(kk,2,2,-i)= eenewtor(kk,2,2,i)
 882         enddo
 883       enddo
 884       if (lprint) then
 885         write (iout,'(a)')
 886      &    "Single-body coefficients of the torsional potentials"
 887         do i=-nloctyp+1,nloctyp-1
 888           write (iout,*) "Type: ",onelet(iloctyp(i))
 889           write (iout,*) "Coefficients of the expansion of B1tor"
 890           do j=1,2
 891             write (iout,'(3hB1(,i1,1h),3f10.5)')
 892      &        j,(bnew1tor(k,j,i),k=1,3)
 893           enddo
 894           write (iout,*) "Coefficients of the expansion of B2tor"
 895           do j=1,2
 896             write (iout,'(3hB2(,i1,1h),3f10.5)')
 897      &        j,(bnew2tor(k,j,i),k=1,3)
 898           enddo
 899           write (iout,*) "Coefficients of the expansion of Ctor"
 900           write (iout,'(3hC11,3f10.5)') (ccnewtor(j,1,i),j=1,3)
 901           write (iout,'(3hC12,3f10.5)') (ccnewtor(j,2,i),j=1,3)
 902           write (iout,*) "Coefficients of the expansion of Dtor"
 903           write (iout,'(3hD11,3f10.5)') (ddnewtor(j,1,i),j=1,3)
 904           write (iout,'(3hD12,3f10.5)') (ddnewtor(j,2,i),j=1,3)
 905           write (iout,*) "Coefficients of the expansion of Etor"
 906           write (iout,'(2hE0,3f10.5)') (e0newtor(j,i),j=1,3)
 907           do j=1,2
 908             do k=1,2
 909               write (iout,'(1hE,2i1,2f10.5)')
 910      &          j,k,(eenewtor(l,j,k,i),l=1,2)
 911             enddo
 912           enddo
 913         enddo
 914       endif
 915       ELSE
 916       do i=-nloctyp+1,nloctyp-1
 917         do ii=1,3
 918           do j=1,2
 919             bnew1tor(ii,j,i)=bnew1(ii,j,i)
 920           enddo
 921         enddo
 922         do ii=1,3
 923           do j=1,2
 924             bnew2tor(ii,j,i)=bnew2(ii,j,i)
 925           enddo
 926         enddo
 927         do ii=1,3
 928           ccnewtor(ii,1,i)=ccnew(ii,1,i)
 929           ccnewtor(ii,2,i)=ccnew(ii,2,i)
 930           ddnewtor(ii,1,i)=ddnew(ii,1,i)
 931           ddnewtor(ii,2,i)=ddnew(ii,2,i)
 932         enddo
 933       enddo
 934       ENDIF !SPLIT_FOURIER_TOR
 935 #else
 936       if (lprint)
 937      &  write (iout,*) "Coefficients of the expansion of Eloc(l1,l2)"
 938       do i=0,nloctyp-1
 939         read (ifourier,*,end=115,err=115)
 940         read (ifourier,*,end=115,err=115) (b(ii,i),ii=1,13)
 941         if (lprint) then
 942         write (iout,*) 'Type ',onelet(iloctyp(i))
 943         write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
 944         endif
 945         if (i.gt.0) then
 946         b(2,-i)= b(2,i)
 947         b(3,-i)= b(3,i)
 948         b(4,-i)=-b(4,i)
 949         b(5,-i)=-b(5,i)
 950         endif
 951 c        B1(1,i)  = b(3)
 952 c        B1(2,i)  = b(5)
 953 c        B1(1,-i) = b(3)
 954 c        B1(2,-i) = -b(5)
 955 c        b1(1,i)=0.0d0
 956 c        b1(2,i)=0.0d0
 957 c        B1tilde(1,i) = b(3)
 958 c        B1tilde(2,i) =-b(5)
 959 c        B1tilde(1,-i) =-b(3)
 960 c        B1tilde(2,-i) =b(5)
 961 c        b1tilde(1,i)=0.0d0
 962 c        b1tilde(2,i)=0.0d0
 963 c        B2(1,i)  = b(2)
 964 c        B2(2,i)  = b(4)
 965 c        B2(1,-i)  =b(2)
 966 c        B2(2,-i)  =-b(4)
 967 cc        B1tilde(1,i) = b(3,i)
 968 cc        B1tilde(2,i) =-b(5,i)
 969 C        B1tilde(1,-i) =-b(3,i)
 970 C        B1tilde(2,-i) =b(5,i)
 971 cc        b1tilde(1,i)=0.0d0
 972 cc        b1tilde(2,i)=0.0d0
 973 cc        B2(1,i)  = b(2,i)
 974 cc        B2(2,i)  = b(4,i)
 975 C        B2(1,-i)  =b(2,i)
 976 C        B2(2,-i)  =-b(4,i)
 977
 978 c        b2(1,i)=0.0d0
 979 c        b2(2,i)=0.0d0
 980         CCold(1,1,i)= b(7,i)
 981         CCold(2,2,i)=-b(7,i)
 982         CCold(2,1,i)= b(9,i)
 983         CCold(1,2,i)= b(9,i)
 984         CCold(1,1,-i)= b(7,i)
 985         CCold(2,2,-i)=-b(7,i)
 986         CCold(2,1,-i)=-b(9,i)
 987         CCold(1,2,-i)=-b(9,i)
 988 c        CC(1,1,i)=0.0d0
 989 c        CC(2,2,i)=0.0d0
 990 c        CC(2,1,i)=0.0d0
 991 c        CC(1,2,i)=0.0d0
 992 c        Ctilde(1,1,i)= CCold(1,1,i)
 993 c        Ctilde(1,2,i)= CCold(1,2,i)
 994 c        Ctilde(2,1,i)=-CCold(2,1,i)
 995 c        Ctilde(2,2,i)=-CCold(2,2,i)
 996
 997 c        Ctilde(1,1,i)=0.0d0
 998 c        Ctilde(1,2,i)=0.0d0
 999 c        Ctilde(2,1,i)=0.0d0
1000 c        Ctilde(2,2,i)=0.0d0
1001         DDold(1,1,i)= b(6,i)
1002         DDold(2,2,i)=-b(6,i)
1003         DDold(2,1,i)= b(8,i)
1004         DDold(1,2,i)= b(8,i)
1005         DDold(1,1,-i)= b(6,i)
1006         DDold(2,2,-i)=-b(6,i)
1007         DDold(2,1,-i)=-b(8,i)
1008         DDold(1,2,-i)=-b(8,i)
1009 c        DD(1,1,i)=0.0d0
1010 c        DD(2,2,i)=0.0d0
1011 c        DD(2,1,i)=0.0d0
1012 c        DD(1,2,i)=0.0d0
1013 c        Dtilde(1,1,i)= DD(1,1,i)
1014 c        Dtilde(1,2,i)= DD(1,2,i)
1015 c        Dtilde(2,1,i)=-DD(2,1,i)
1016 c        Dtilde(2,2,i)=-DD(2,2,i)
1017
1018 c        Dtilde(1,1,i)=0.0d0
1019 c        Dtilde(1,2,i)=0.0d0
1020 c        Dtilde(2,1,i)=0.0d0
1021 c        Dtilde(2,2,i)=0.0d0
1022         EEold(1,1,i)= b(10,i)+b(11,i)
1023         EEold(2,2,i)=-b(10,i)+b(11,i)
1024         EEold(2,1,i)= b(12,i)-b(13,i)
1025         EEold(1,2,i)= b(12,i)+b(13,i)
1026         EEold(1,1,-i)= b(10,i)+b(11,i)
1027         EEold(2,2,-i)=-b(10,i)+b(11,i)
1028         EEold(2,1,-i)=-b(12,i)+b(13,i)
1029         EEold(1,2,-i)=-b(12,i)-b(13,i)
1030         write(iout,*) "TU DOCHODZE"
1031         print *,"JESTEM"
1032 c        ee(1,1,i)=1.0d0
1033 c        ee(2,2,i)=1.0d0
1034 c        ee(2,1,i)=0.0d0
1035 c        ee(1,2,i)=0.0d0
1036 c        ee(2,1,i)=ee(1,2,i)
1037       enddo
1038       if (lprint) then
1039       write (iout,*)
1040       write (iout,*)
1041      &"Coefficients of the cumulants (independent of valence angles)"
1042       do i=-nloctyp+1,nloctyp-1
1043         write (iout,*) 'Type ',onelet(iloctyp(i))
1044         write (iout,*) 'B1'
1045         write(iout,'(2f10.5)') B(3,i),B(5,i)
1046         write (iout,*) 'B2'
1047         write(iout,'(2f10.5)') B(2,i),B(4,i)
1048         write (iout,*) 'CC'
1049         do j=1,2
1050           write (iout,'(2f10.5)') CCold(j,1,i),CCold(j,2,i)
1051         enddo
1052         write(iout,*) 'DD'
1053         do j=1,2
1054           write (iout,'(2f10.5)') DDold(j,1,i),DDold(j,2,i)
1055         enddo
1056         write(iout,*) 'EE'
1057         do j=1,2
1058           write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i)
1059         enddo
1060       enddo
1061       endif
1062 #endif
1063 C      write (iout,*) 'KURWAKURWA'
1064 #ifdef CRYST_TOR
1065 C
1066 C Read torsional parameters in old format
1067 C
1068       read (itorp,*,end=113,err=113) ntortyp,nterm_old
1069       write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
1070       read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1071       do i=1,ntortyp
1072         do j=1,ntortyp
1073           read (itorp,'(a)',end=113,err=113)
1074           do k=1,nterm_old
1075             read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
1076           enddo
1077         enddo
1078       enddo
1079       close (itorp)
1080       if (lprint) then
1081         write (iout,'(/a/)') 'Torsional constants:'
1082         do i=1,ntortyp
1083           do j=1,ntortyp
1084             write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
1085             write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
1086           enddo
1087         enddo
1088       endif
1089 #else
1090 C
1091 C Read torsional parameters
1092 C
1093       IF (TOR_MODE.eq.0) THEN
1094
1095       read (itorp,*,end=113,err=113) ntortyp
1096       read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1097       write (iout,*) 'ntortyp',ntortyp
1098       do i=1,ntyp1
1099         itype2loc(i)=itortyp(i)
1100       enddo
1101       do i=1,ntyp1
1102         itype2loc(-i)=-itype2loc(i)
1103       enddo
1104       do iblock=1,2
1105       do i=-ntyp,-1
1106        itortyp(i)=-itortyp(-i)
1107       enddo
1108 c      write (iout,*) 'ntortyp',ntortyp
1109       do i=0,ntortyp-1
1110         do j=-ntortyp+1,ntortyp-1
1111           read (itorp,*,end=113,err=113) nterm(i,j,iblock),
1112      &          nlor(i,j,iblock)
1113           nterm(-i,-j,iblock)=nterm(i,j,iblock)
1114           nlor(-i,-j,iblock)=nlor(i,j,iblock)
1115           v0ij=0.0d0
1116           si=-1.0d0
1117           do k=1,nterm(i,j,iblock)
1118             read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),
1119      &      v2(k,i,j,iblock)
1120             v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
1121             v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
1122             v0ij=v0ij+si*v1(k,i,j,iblock)
1123             si=-si
1124          enddo
1125           do k=1,nlor(i,j,iblock)
1126             read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),
1127      &        vlor2(k,i,j),vlor3(k,i,j)
1128             v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
1129           enddo
1130           v0(i,j,iblock)=v0ij
1131           v0(-i,-j,iblock)=v0ij
1132         enddo
1133       enddo
1134       enddo
1135       close (itorp)
1136       if (lprint) then
1137         write (iout,'(/a/)') 'Torsional constants:'
1138         do iblock=1,2
1139         do i=1,ntortyp
1140           do j=1,ntortyp
1141             write (iout,*) 'ityp',i,' jtyp',j
1142             write (iout,*) 'Fourier constants'
1143             do k=1,nterm(i,j,iblock)
1144               write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
1145      &        v2(k,i,j,iblock)
1146             enddo
1147             write (iout,*) 'Lorenz constants'
1148             do k=1,nlor(i,j,iblock)
1149               write (iout,'(3(1pe15.5))')
1150      &         vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
1151             enddo
1152           enddo
1153         enddo
1154         enddo
1155       endif
1156 C
1157 C 6/23/01 Read parameters for double torsionals
1158 C
1159       do iblock=1,2
1160       do i=0,ntortyp-1
1161         do j=-ntortyp+1,ntortyp-1
1162           do k=-ntortyp+1,ntortyp-1
1163             read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
1164 c              write (iout,*) "OK onelett",
1165 c     &         i,j,k,t1,t2,t3
1166
1167             if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
1168      &        .or. t3.ne.toronelet(k)) then
1169               write (iout,*) "Error in double torsional parameter file",
1170      &         i,j,k,t1,t2,t3
1171 #ifdef MPI
1172               call MPI_Finalize(Ierror)
1173 #endif
1174                stop "Error in double torsional parameter file"
1175             endif
1176           read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),
1177      &         ntermd_2(i,j,k,iblock)
1178             ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
1179             ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
1180             read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,
1181      &         ntermd_1(i,j,k,iblock))
1182             read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,
1183      &         ntermd_1(i,j,k,iblock))
1184             read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,
1185      &         ntermd_1(i,j,k,iblock))
1186             read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,
1187      &         ntermd_1(i,j,k,iblock))
1188 C Martix of D parameters for one dimesional foureir series
1189             do l=1,ntermd_1(i,j,k,iblock)
1190              v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
1191              v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
1192              v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
1193              v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
1194 c            write(iout,*) "whcodze" ,
1195 c     & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
1196             enddo
1197             read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),
1198      &         v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
1199      &         v2s(m,l,i,j,k,iblock),
1200      &         m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
1201 C Martix of D parameters for two dimesional fourier series
1202             do l=1,ntermd_2(i,j,k,iblock)
1203              do m=1,l-1
1204              v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
1205              v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
1206              v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
1207              v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
1208              enddo!m
1209             enddo!l
1210           enddo!k
1211         enddo!j
1212       enddo!i
1213       enddo!iblock
1214       if (lprint) then
1215       write (iout,*)
1216       write (iout,*) 'Constants for double torsionals'
1217       do iblock=1,2
1218       do i=0,ntortyp-1
1219         do j=-ntortyp+1,ntortyp-1
1220           do k=-ntortyp+1,ntortyp-1
1221             write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
1222      &        ' nsingle',ntermd_1(i,j,k,iblock),
1223      &        ' ndouble',ntermd_2(i,j,k,iblock)
1224             write (iout,*)
1225             write (iout,*) 'Single angles:'
1226             do l=1,ntermd_1(i,j,k,iblock)
1227               write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
1228      &           v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
1229      &           v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),
1230      &           v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
1231             enddo
1232             write (iout,*)
1233             write (iout,*) 'Pairs of angles:'
1234             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1235             do l=1,ntermd_2(i,j,k,iblock)
1236               write (iout,'(i5,20f10.5)')
1237      &         l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
1238             enddo
1239             write (iout,*)
1240            write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1241             do l=1,ntermd_2(i,j,k,iblock)
1242               write (iout,'(i5,20f10.5)')
1243      &         l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),
1244      &         (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
1245             enddo
1246             write (iout,*)
1247           enddo
1248         enddo
1249       enddo
1250       enddo
1251       endif
1252
1253       ELSE IF (TOR_MODE.eq.1) THEN
1254
1255 C read valence-torsional parameters
1256       read (itorp,*,end=113,err=113) ntortyp
1257       nkcctyp=ntortyp
1258       write (iout,*) "Valence-torsional parameters read in ntortyp",
1259      &   ntortyp
1260       read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1261       write (iout,*) "itortyp_kcc",(itortyp(i),i=1,ntyp)
1262       do i=-ntyp,-1
1263         itortyp(i)=-itortyp(-i)
1264       enddo
1265       do i=-ntortyp+1,ntortyp-1
1266         do j=-ntortyp+1,ntortyp-1
1267 C first we read the cos and sin gamma parameters
1268           read (itorp,'(13x,a)',end=113,err=113) string
1269           write (iout,*) i,j,string
1270           read (itorp,*,end=113,err=113)
1271      &    nterm_kcc(j,i),nterm_kcc_Tb(j,i)
1272 C           read (itorkcc,*,end=121,err=121) nterm_kcc_Tb(j,i)
1273           do k=1,nterm_kcc(j,i)
1274             do l=1,nterm_kcc_Tb(j,i)
1275               do ll=1,nterm_kcc_Tb(j,i)
1276               read (itorp,*,end=113,err=113) ii,jj,kk,
1277      &          v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
1278               enddo
1279             enddo
1280           enddo
1281         enddo
1282       enddo
1283
1284       ELSE
1285 c AL 4/8/16: Calculate coefficients from one-body parameters
1286       ntortyp=nloctyp
1287       do i=-ntyp1,ntyp1
1288         itortyp(i)=itype2loc(i)
1289 c        itortyp(i)=0
1290       enddo
1291       write (iout,*)
1292      &"Val-tor parameters calculated from cumulant coefficients ntortyp"
1293      & ,ntortyp
1294       do i=-ntortyp+1,ntortyp-1
1295         do j=-ntortyp+1,ntortyp-1
1296           nterm_kcc(j,i)=2
1297           nterm_kcc_Tb(j,i)=3
1298           do k=1,nterm_kcc_Tb(j,i)
1299             do l=1,nterm_kcc_Tb(j,i)
1300               v1_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,1,j)
1301      &                         +bnew1tor(k,2,i)*bnew2tor(l,2,j)
1302               v2_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,2,j)
1303      &                         +bnew1tor(k,2,i)*bnew2tor(l,1,j)
1304             enddo
1305           enddo
1306           do k=1,nterm_kcc_Tb(j,i)
1307             do l=1,nterm_kcc_Tb(j,i)
1308 c              v1_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,1,i)*ddnewtor(l,1,j)
1309 c     &                         -ccnewtor(k,2,i)*ddnewtor(l,2,j))
1310 c              v2_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,2,i)*ddnewtor(l,1,j)
1311 c     &                         +ccnewtor(k,1,i)*ddnewtor(l,2,j))
1312               v1_kcc(k,l,2,i,j)=-(ccnewtor(k,1,i)*ddnewtor(l,1,j)
1313      &                         -ccnewtor(k,2,i)*ddnewtor(l,2,j))
1314               v2_kcc(k,l,2,i,j)=-(ccnewtor(k,2,i)*ddnewtor(l,1,j)
1315      &                         +ccnewtor(k,1,i)*ddnewtor(l,2,j))
1316             enddo
1317           enddo
1318 cf(theta,gamma)=-(b21(theta)*b11(theta)+b12(theta)*b22(theta))*cos(gamma)-(b22(theta)*b11(theta)+b21(theta)*b12(theta))*sin(gamma)+(c11(theta)*d11(theta)-c12(theta)*d12(theta))*cos(2*gamma)+(c12(theta)*d11(theta)+c11(theta)*d12(theta))*sin(2*gamma)
1319         enddo
1320       enddo
1321
1322       ENDIF ! TOR_MODE
1323
1324       if (tor_mode.gt.0 .and. lprint) then
1325 c Print valence-torsional parameters
1326         write (iout,'(a)')
1327      &    "Parameters of the valence-torsional potentials"
1328         do i=-ntortyp+1,ntortyp-1
1329         do j=-ntortyp+1,ntortyp-1
1330         write (iout,'(3a)') "Type ",toronelet(i),toronelet(j)
1331         write (iout,'(3a5,2a15)') "itor","ival","jval","v_kcc","v2_kcc"
1332         do k=1,nterm_kcc(j,i)
1333           do l=1,nterm_kcc_Tb(j,i)
1334             do ll=1,nterm_kcc_Tb(j,i)
1335                write (iout,'(3i5,2f15.4)')
1336      &            k,l-1,ll-1,v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
1337             enddo
1338           enddo
1339         enddo
1340         enddo
1341         enddo
1342       endif
1343
1344 #endif
1345 C Read of Side-chain backbone correlation parameters
1346 C Modified 11 May 2012 by Adasko
1347 CCC
1348 C
1349       read (isccor,*,end=119,err=119) nsccortyp
1350       read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
1351       do i=-ntyp,-1
1352         isccortyp(i)=-isccortyp(-i)
1353       enddo
1354       iscprol=isccortyp(20)
1355 c      write (iout,*) 'ntortyp',ntortyp
1356       maxinter=3
1357 cc maxinter is maximum interaction sites
1358       do l=1,maxinter
1359       do i=1,nsccortyp
1360         do j=1,nsccortyp
1361           read (isccor,*,end=119,err=119)
1362      &nterm_sccor(i,j),nlor_sccor(i,j)
1363 c          write (iout,*) nterm_sccor(i,j)
1364           v0ijsccor=0.0d0
1365           v0ijsccor1=0.0d0
1366           v0ijsccor2=0.0d0
1367           v0ijsccor3=0.0d0
1368           si=-1.0d0
1369           nterm_sccor(-i,j)=nterm_sccor(i,j)
1370           nterm_sccor(-i,-j)=nterm_sccor(i,j)
1371           nterm_sccor(i,-j)=nterm_sccor(i,j)
1372 c          write (iout,*) nterm_sccor(i,j),nterm_sccor(-i,j),
1373 c     &    nterm_sccor(-i,-j),nterm_sccor(i,-j)
1374           do k=1,nterm_sccor(i,j)
1375             read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
1376      &    ,v2sccor(k,l,i,j)
1377             if (j.eq.iscprol) then
1378              if (i.eq.isccortyp(10)) then
1379              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1380              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1381              else
1382              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
1383      &                        +v2sccor(k,l,i,j)*dsqrt(0.75d0)
1384              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
1385      &                        +v1sccor(k,l,i,j)*dsqrt(0.75d0)
1386              v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1387              v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1388              v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1389              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1390              endif
1391             else
1392              if (i.eq.isccortyp(10)) then
1393              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1394              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1395              else
1396                if (j.eq.isccortyp(10)) then
1397              v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
1398              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
1399                else
1400              v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
1401              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1402              v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1403              v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1404              v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1405              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1406                 endif
1407                endif
1408             endif
1409             v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
1410             v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
1411             v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
1412             v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
1413             si=-si
1414            enddo
1415           do k=1,nlor_sccor(i,j)
1416             read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
1417      &        vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1418             v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
1419      &(1+vlor3sccor(k,i,j)**2)
1420           enddo
1421           v0sccor(l,i,j)=v0ijsccor
1422           v0sccor(l,-i,j)=v0ijsccor1
1423           v0sccor(l,i,-j)=v0ijsccor2
1424           v0sccor(l,-i,-j)=v0ijsccor3
1425           enddo
1426         enddo
1427       enddo
1428       close (isccor)
1429       if (lprint) then
1430         write (iout,'(/a/)') 'Torsional constants of SCCORR:'
1431         do l=1,maxinter
1432         do i=1,nsccortyp
1433           do j=1,nsccortyp
1434             write (iout,*) 'ityp',i,' jtyp',j
1435             write (iout,*) 'Fourier constants'
1436             do k=1,nterm_sccor(i,j)
1437               write (iout,'(2(1pe15.5))')
1438      & v1sccor(k,l,i,j),v2sccor(k,l,i,j)
1439             enddo
1440             write (iout,*) 'Lorenz constants'
1441             do k=1,nlor_sccor(i,j)
1442               write (iout,'(3(1pe15.5))')
1443      &         vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1444             enddo
1445           enddo
1446         enddo
1447         enddo
1448       endif
1449 C
1450 C Read electrostatic-interaction parameters
1451 C
1452       if (lprint) then
1453         write (iout,'(/a)') 'Electrostatic interaction constants:'
1454         write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
1455      &            'IT','JT','APP','BPP','AEL6','AEL3'
1456       endif
1457       read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
1458       read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
1459       read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
1460       read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
1461       close (ielep)
1462       do i=1,2
1463         do j=1,2
1464         rri=rpp(i,j)**6
1465         app (i,j)=epp(i,j)*rri*rri
1466         bpp (i,j)=-2.0D0*epp(i,j)*rri
1467         ael6(i,j)=elpp6(i,j)*4.2D0**6
1468         ael3(i,j)=elpp3(i,j)*4.2D0**3
1469         lprint=.true.
1470         if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
1471      &                    ael6(i,j),ael3(i,j)
1472         enddo
1473       enddo
1474 C
1475 C Read side-chain interaction parameters.
1476 C
1477       read (isidep,*) ipot,expon
1478       if (ipot.lt.1 .or. ipot.gt.6) then
1479         write (iout,'(2a)') 'Error while reading SC interaction',
1480      &               'potential file - unknown potential type.'
1481         stop
1482       endif
1483       expon2=expon/2
1484       write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
1485      & ', exponents are ',expon,2*expon
1486       goto (10,20,30,30,40,50,60) ipot
1487 C----------------------- LJ potential ---------------------------------
1488    10 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
1489       if (lprint) then
1490         write (iout,'(/a/)') 'Parameters of the LJ potential:'
1491         write (iout,'(a/)') 'The epsilon array:'
1492         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1493         write (iout,'(/a)') 'One-body parameters:'
1494         write (iout,'(a,4x,a)') 'residue','sigma'
1495         write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1496       endif
1497       goto 60
1498 C----------------------- LJK potential --------------------------------
1499    20 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
1500      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1501       if (lprint) then
1502         write (iout,'(/a/)') 'Parameters of the LJK potential:'
1503         write (iout,'(a/)') 'The epsilon array:'
1504         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1505         write (iout,'(/a)') 'One-body parameters:'
1506         write (iout,'(a,4x,2a)') 'residue','   sigma  ','    r0    '
1507         write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
1508      &        i=1,ntyp)
1509       endif
1510       goto 60
1511 C---------------------- GB or BP potential -----------------------------
1512    30 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
1513      &  (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip0(i),i=1,ntyp),
1514      &  (alp(i),i=1,ntyp)
1515 C For the GB potential convert sigma'**2 into chi'
1516       do i=1,ntyp
1517         chip(i)=(chip0(i)-1.0D0)/(chip0(i)+1.0D0)
1518       enddo
1519       if (lprint) then
1520         write (iout,'(/a/)') 'Parameters of the BP/GB potential:'
1521         write (iout,'(a/)') 'The epsilon array:'
1522         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1523         write (iout,'(/a)') 'One-body parameters:'
1524         write (iout,'(a,4x,5a)') 'residue','   sigma  ','s||/s_|_^2',
1525      &       '   chip0  ','    chip  ','    alph  '
1526         write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),sigii(i),
1527      &                     chip0(i),chip(i),alp(i),i=1,ntyp)
1528       endif
1529       goto 60
1530 C--------------------- GBV potential -----------------------------------
1531 c   40 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
1532 c     &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1533 c     &  (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1534    40 do i=1,ntyp
1535        read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
1536       enddo
1537       read (isidep,*,end=117,err=117)(sigma0(i),i=1,ntyp)
1538       read (isidep,*,end=117,err=117)(rr0(i),i=1,ntyp)
1539       read (isidep,*,end=117,err=117)(sigii(i),i=1,ntyp)
1540       read (isidep,*,end=117,err=117)(chip0(i),i=1,ntyp)
1541       read (isidep,*,end=117,err=117)(alp(i),i=1,ntyp)
1542 c      do i=1,ntyp
1543 c       read (isidep,*,end=1161,err=1161)(epslip(i,j),j=i,ntyp)
1544 c       if (lprint) write(iout,*) "epslip", i, (epslip(i,j),j=i,ntyp)
1545 c      enddo
1546       do i=1,ntyp
1547         chip(i)=(chip0(i)-1.0D0)/(chip0(i)+1.0D0)
1548       enddo
1549       if (lprint) then
1550         write (iout,'(/a/)') 'Parameters of the GBV potential:'
1551         write (iout,'(a/)') 'The epsilon array:'
1552         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1553         write (iout,'(/a)') 'One-body parameters:'
1554         write (iout,'(a,4x,6a)') 'residue','   sigma  ','    r0    ',
1555      &      's||/s_|_^2','   chip0  ','    chip  ','    alph  '
1556         write (iout,'(a3,6x,6f10.5)') (restyp(i),sigma0(i),rr0(i),
1557      &           sigii(i),chip0(i),chip(i),alp(i),i=1,ntyp)
1558       endif
1559       goto 60
1560 C--------------------- Momo potential -----------------------------------
1561    50 continue
1562       read (isidep,*) (icharge(i),i=1,ntyp)
1563 c      write (2,*) "icharge",(icharge(i),i=1,ntyp)
1564       do i=1,ntyp
1565        do j=1,i
1566 c!        write (*,*) "Im in ", i, " ", j
1567         read(isidep,*)
1568      &  eps(i,j),sigma(i,j),chi(i,j),chi(j,i),chipp(i,j),chipp(j,i),
1569      &  (alphasur(k,i,j),k=1,4),sigmap1(i,j),sigmap2(i,j),
1570      &  chis(i,j),chis(j,i),
1571      &  nstate(i,j),(wstate(k,i,j),k=1,4),
1572      &  dhead(1,1,i,j),
1573      &  dhead(1,2,i,j),
1574      &  dhead(2,1,i,j),
1575      &  dhead(2,2,i,j),
1576      &  dtail(1,i,j),dtail(2,i,j),
1577      &  epshead(i,j),sig0head(i,j),
1578      &  rborn(i,j),rborn(j,i),
1579      &  (wqdip(k,i,j),k=1,2),wquad(i,j),
1580      &  alphapol(i,j),alphapol(j,i),
1581      &  (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(i,j),epsintab(i,j)
1582        END DO
1583       END DO
1584 c!      write (*,*) "nstate gly-gly", nstate(10,10)
1585 c! THIS LOOP FILLS PARAMETERS FOR PAIRS OF AA's NOT EXPLICITLY
1586 c! DEFINED IN SCPARM.MOMO. IT DOES SO BY TAKING THEM FROM SYMMETRIC
1587 c! PAIR, FOR EG. IF ARG-HIS IS BLANK, IT WILL TAKE PARAMETERS
1588 c! FROM HIS-ARG.
1589 c!
1590 c! DISABLE IT AT >>YOUR OWN PERIL<<
1591 c!
1592       DO i = 1, ntyp
1593        DO j = i+1, ntyp
1594         eps(i,j) = eps(j,i)
1595         sigma(i,j) = sigma(j,i)
1596         nstate(i,j) = nstate(j,i)
1597         sigmap1(i,j) = sigmap1(j,i)
1598         sigmap2(i,j) = sigmap2(j,i)
1599         sigiso1(i,j) = sigiso1(j,i)
1600         sigiso2(i,j) = sigiso2(j,i)
1601
1602         DO k = 1, 4
1603          alphasur(k,i,j) = alphasur(k,j,i)
1604          wstate(k,i,j) = wstate(k,j,i)
1605          alphiso(k,i,j) = alphiso(k,j,i)
1606         END DO
1607
1608         dhead(2,1,i,j) = dhead(1,1,j,i)
1609         dhead(2,2,i,j) = dhead(1,2,j,i)
1610         dhead(1,1,i,j) = dhead(2,1,j,i)
1611         dhead(1,2,i,j) = dhead(2,2,j,i)
1612         dtail(1,i,j) = dtail(1,j,i)
1613         dtail(2,i,j) = dtail(2,j,i)
1614 c!        DO k = 1, 2
1615 c!         DO l = 1, 2
1616 c!         write (*,*) "dhead(k,l,j,i) = ", dhead(k,l,j,i)
1617 c!         write (*,*) "dhead(k,l,i,j) = ", dhead(k,l,i,j)
1618 c!          dhead(l,k,i,j) = dhead(k,l,j,i)
1619 c!         END DO
1620 c!        END DO
1621
1622         epshead(i,j) = epshead(j,i)
1623         sig0head(i,j) = sig0head(j,i)
1624
1625         DO k = 1, 2
1626          wqdip(k,i,j) = wqdip(k,j,i)
1627         END DO
1628
1629         wquad(i,j) = wquad(j,i)
1630         epsintab(i,j) = epsintab(j,i)
1631
1632        END DO
1633       END DO
1634
1635       if (.not.lprint) goto 70
1636       write (iout,'(a)')
1637      & "Parameters of the new physics-based SC-SC interaction potential"
1638       write (iout,'(/7a)') 'Residues','     epsGB','       rGB',
1639      &  '    chi1GB','    chi2GB','   chip1GB','   chip2GB'
1640       do i=1,ntyp
1641         do j=1,i
1642           write (iout,'(2(a3,1x),1pe10.3,5(0pf10.3))')
1643      &      restyp(i),restyp(j),eps(i,j),sigma(i,j),chi(i,j),chi(j,i),
1644      &       chipp(i,j),chipp(j,i)
1645         enddo
1646       enddo
1647       write (iout,'(/9a)') 'Residues',' alphasur1',' alphasur2',
1648      &  ' alphasur3',' alphasur4','   sigmap1','   sigmap2',
1649      &  '     chis1','     chis2'
1650       do i=1,ntyp
1651         do j=1,i
1652           write (iout,'(2(a3,1x),8(0pf10.3))')
1653      &      restyp(i),restyp(j),(alphasur(k,i,j),k=1,4),
1654      &       sigmap1(i,j),sigmap2(j,i),chis(i,j),chis(j,i)
1655         enddo
1656       enddo
1657       write (iout,'(/14a)') 'Residues',' nst ','  wst1',
1658      &  '    wst2','    wst3','    wst4','   dh11','   dh21',
1659      &  '   dh12','   dh22','    dt1','    dt2','    epsh1',
1660      &  '   sigh'
1661       do i=1,ntyp
1662         do j=1,i
1663           write (iout,'(2(a3,1x),i3,4f8.3,6f7.2,f9.5,f7.2)')
1664      &      restyp(i),restyp(j),nstate(i,j),(wstate(k,i,j),k=1,4),
1665      &       ((dhead(l,k,i,j),l=1,2),k=1,2),dtail(1,i,j),dtail(2,i,j),
1666      &       epshead(i,j),sig0head(i,j)
1667         enddo
1668       enddo
1669       write (iout,'(/12a)') 'Residues',' ch1',' ch2',
1670      &  '    rborn1','    rborn2','    wqdip1','    wqdip2',
1671      &  '     wquad'
1672       do i=1,ntyp
1673         do j=1,i
1674           write (iout,'(2(a3,1x),2i4,5f10.3)')
1675      &      restyp(i),restyp(j),icharge(i),icharge(j),
1676      &      rborn(i,j),rborn(j,i),(wqdip(k,i,j),k=1,2),wquad(i,j)
1677         enddo
1678       enddo
1679       write (iout,'(/12a)') 'Residues',
1680      &  '  alphpol1',
1681      &  '  alphpol2','  alphiso1','   alpiso2',
1682      &  '   alpiso3','   alpiso4','   sigiso1','   sigiso2',
1683      &  '     epsin'
1684       do i=1,ntyp
1685         do j=1,i
1686           write (iout,'(2(a3,1x),11f10.3)')
1687      &      restyp(i),restyp(j),alphapol(i,j),alphapol(j,i),
1688      &      (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(j,i),
1689      &      epsintab(i,j)
1690         enddo
1691       enddo
1692       goto 70
1693 C For the GB potential convert sigma'**2 into chi'
1694       DO i=1,ntyp
1695        chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1696       END DO
1697       IF (lprint) THEN
1698        write (iout,'(/a/)') 'Parameters of the BP potential:'
1699        write (iout,'(a/)') 'The epsilon array:'
1700        CALL printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1701        write (iout,'(/a)') 'One-body parameters:'
1702        write (iout,'(a,4x,4a)') 'residue','   sigma  ','s||/s_|_^2',
1703      &       '    chip  ','    alph  '
1704        write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1705      &                     chip(i),alp(i),i=1,ntyp)
1706       END IF
1707    60 CONTINUE
1708       close (isidep)
1709 C-----------------------------------------------------------------------
1710 C Calculate the "working" parameters of SC interactions.
1711       do i=2,ntyp
1712         do j=1,i-1
1713         eps(i,j)=eps(j,i)
1714           do k=1,4
1715             alphasur(k,j,i)=alphasur(k,i,j)
1716             alphiso(k,j,i)=alphiso(k,i,j)
1717             wstate(k,j,i)=wstate(k,i,j)
1718           enddo
1719           do k=1,2
1720             wqdip(k,j,i)=wqdip(k,i,j)
1721           enddo
1722           do k=1,2
1723             do l=1,2
1724               dhead(l,k,j,i)=dhead(l,k,i,j)
1725             enddo
1726           enddo
1727         enddo
1728       enddo
1729       IF (ipot.LT.6) THEN
1730        do i=1,ntyp
1731         do j=i,ntyp
1732          sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1733          sigma(j,i)=sigma(i,j)
1734          rs0(i,j)=dwa16*sigma(i,j)
1735          rs0(j,i)=rs0(i,j)
1736         enddo
1737        enddo
1738       END IF
1739
1740    70 continue
1741       if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1742      & 'Working parameters of the SC interactions:',
1743      & '     a    ','     b    ','   augm   ','  sigma ','   r0   ',
1744      & '  chi1   ','   chi2   '
1745       do i=1,ntyp
1746         do j=i,ntyp
1747           epsij=eps(i,j)
1748           if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4 .or. ipot.eq.6 )
1749      &    THEN
1750             rrij=sigma(i,j)
1751           else
1752             rrij=rr0(i)+rr0(j)
1753           endif
1754           r0(i,j)=rrij
1755           r0(j,i)=rrij
1756           rrij=rrij**expon
1757           epsij=eps(i,j)
1758           sigeps=dsign(1.0D0,epsij)
1759           epsij=dabs(epsij)
1760           aa(i,j)=epsij*rrij*rrij
1761           bb(i,j)=-sigeps*epsij*rrij
1762           aa(j,i)=aa(i,j)
1763           bb(j,i)=bb(i,j)
1764           IF ((ipot.gt.2).AND.(ipot.LT.6)) THEN
1765             sigt1sq=sigma0(i)**2
1766             sigt2sq=sigma0(j)**2
1767             sigii1=sigii(i)
1768             sigii2=sigii(j)
1769             ratsig1=sigt2sq/sigt1sq
1770             ratsig2=1.0D0/ratsig1
1771             chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1772             if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1773             rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1774           else
1775             rsum_max=sigma(i,j)
1776           endif
1777 c         if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1778             sigmaii(i,j)=rsum_max
1779             sigmaii(j,i)=rsum_max
1780 c         else
1781 c           sigmaii(i,j)=r0(i,j)
1782 c           sigmaii(j,i)=r0(i,j)
1783 c         endif
1784 cd        write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1785           if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1786             r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1787             augm(i,j)=epsij*r_augm**(2*expon)
1788 c           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1789             augm(j,i)=augm(i,j)
1790           else
1791             augm(i,j)=0.0D0
1792             augm(j,i)=0.0D0
1793           endif
1794           if (lprint) then
1795             if (ipot.lt.6) then
1796             write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1797      &      restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1798      &      sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1799             else
1800             write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3),2i3,10f8.4,
1801      &       i3,40f10.4)')
1802      &      restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1803      &      sigma(i,j),r0(i,j),chi(i,j),chi(j,i),
1804      &      icharge(i),icharge(j),chipp(i,j),chipp(j,i),
1805      &     (alphasur(k,i,j),k=1,4),sigmap1(i,j),sigmap2(j,i),
1806      &     chis(i,j),chis(j,i),
1807      &     nstate(i,j),(wstate(k,i,j),k=1,4),
1808      &     ((dhead(l,k,i,j),l=1,2),k=1,2),dtail(1,i,j),dtail(2,i,j),
1809      &     epshead(i,j),sig0head(i,j),
1810      &     rborn(i,j),(wqdip(k,i,j),k=1,2),wquad(i,j),
1811      &     alphapol(i,j),alphapol(j,i),
1812      &     (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(i,j)
1813
1814             endif
1815                 endif
1816         enddo
1817       enddo
1818 C
1819 C Define the SC-p interaction constants
1820 C
1821 #ifdef OLDSCP
1822       do i=1,20
1823 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1824 C helix formation)
1825 c       aad(i,1)=0.3D0*4.0D0**12
1826 C Following line for constants currently implemented
1827 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1828         aad(i,1)=1.5D0*4.0D0**12
1829 c       aad(i,1)=0.17D0*5.6D0**12
1830         aad(i,2)=aad(i,1)
1831 C "Soft" SC-p repulsion
1832         bad(i,1)=0.0D0
1833 C Following line for constants currently implemented
1834 c       aad(i,1)=0.3D0*4.0D0**6
1835 C "Hard" SC-p repulsion
1836         bad(i,1)=3.0D0*4.0D0**6
1837 c       bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1838         bad(i,2)=bad(i,1)
1839 c       aad(i,1)=0.0D0
1840 c       aad(i,2)=0.0D0
1841 c       bad(i,1)=1228.8D0
1842 c       bad(i,2)=1228.8D0
1843       enddo
1844 #else
1845 C
1846 C 8/9/01 Read the SC-p interaction constants from file
1847 C
1848       do i=1,ntyp
1849         read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1850       enddo
1851       do i=1,ntyp
1852         aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1853         aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1854         bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1855         bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1856       enddo
1857
1858       if (lprint) then
1859         write (iout,*) "Parameters of SC-p interactions:"
1860         do i=1,20
1861           write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1862      &     eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1863         enddo
1864       endif
1865 #endif
1866 C
1867 C Define the constants of the disulfide bridge
1868 C
1869 C      ebr=-12.0D0
1870 c
1871 c Old arbitrary potential - commented out.
1872 c
1873 c      dbr= 4.20D0
1874 c      fbr= 3.30D0
1875 c
1876 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1877 c energy surface of diethyl disulfide.
1878 c A. Liwo and U. Kozlowska, 11/24/03
1879 c
1880 C      D0CM = 3.78d0
1881 C      AKCM = 15.1d0
1882 C      AKTH = 11.0d0
1883 C      AKCT = 12.0d0
1884 C      V1SS =-1.08d0
1885 C      V2SS = 7.61d0
1886 C      V3SS = 13.7d0
1887       write (iout,*) dyn_ss,'dyndyn'
1888       if (dyn_ss) then
1889         ss_depth=ebr/wsc-0.25*eps(1,1)
1890 C        write(iout,*) akcm,whpb,wsc,'KURWA'
1891         Ht=Ht/wsc-0.25*eps(1,1)
1892
1893         akcm=akcm*whpb/wsc
1894         akth=akth*whpb/wsc
1895         akct=akct*whpb/wsc
1896         v1ss=v1ss*whpb/wsc
1897         v2ss=v2ss*whpb/wsc
1898         v3ss=v3ss*whpb/wsc
1899       else
1900         ss_depth=ebr/whpb-0.25*eps(1,1)*wsc/whpb
1901       endif
1902
1903 C      if (lprint) then
1904       write (iout,'(/a)') "Disulfide bridge parameters:"
1905       write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1906       write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1907       write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1908       write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1909      & ' v3ss:',v3ss
1910 C      endif
1911       if (shield_mode.gt.0) then
1912       pi=3.141592d0
1913 C VSolvSphere the volume of solving sphere
1914 C      print *,pi,"pi"
1915 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
1916 C there will be no distinction between proline peptide group and normal peptide
1917 C group in case of shielding parameters
1918       VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
1919       VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
1920       write (iout,*) VSolvSphere,VSolvSphere_div
1921 C long axis of side chain
1922       do i=1,ntyp
1923       long_r_sidechain(i)=vbldsc0(1,i)
1924       short_r_sidechain(i)=sigma0(i)
1925       enddo
1926       buff_shield=1.0d0
1927       endif
1928       return
1929   111 write (iout,*) "Error reading bending energy parameters."
1930       goto 999
1931   112 write (iout,*) "Error reading rotamer energy parameters."
1932       goto 999
1933   113 write (iout,*) "Error reading torsional energy parameters."
1934       goto 999
1935   114 write (iout,*) "Error reading double torsional energy parameters."
1936       goto 999
1937   115 write (iout,*)
1938      &  "Error reading cumulant (multibody energy) parameters."
1939       goto 999
1940   116 write (iout,*) "Error reading electrostatic energy parameters."
1941       goto 999
1942  1161 write (iout,*) "Error reading lipid energy parameters."
1943       goto 999
1944   117 write (iout,*) "Error reading side chain interaction parameters."
1945       goto 999
1946   118 write (iout,*) "Error reading SCp interaction parameters."
1947       goto 999
1948   119 write (iout,*) "Error reading SCCOR parameters"
1949       goto 999
1950   121 write (iout,*) "Error reading bond parameters"
1951   999 continue
1952 #ifdef MPI
1953       call MPI_Finalize(Ierror)
1954 #endif
1955       stop
1956       return
1957       end