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