source/maxlik/src_MD_T_maxlik-NEWCORR.safe/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),1.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 #ifdef NEWCORR
 679       read (ifourier,*,end=115,err=115) (itype2loc(i),i=1,ntyp)
 680       read (ifourier,*,end=115,err=115) (iloctyp(i),i=0,nloctyp-1)
 681       itype2loc(ntyp1)=nloctyp
 682       iloctyp(nloctyp)=ntyp1
 683       do i=1,ntyp1
 684         itype2loc(-i)=-itype2loc(i)
 685       enddo
 686 #else
 687       do i=1,ntyp1
 688         itype2loc(i)=itortyp(i)
 689       enddo
 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 #else
 809       if (lprint)
 810      &  write (iout,*) "Coefficients of the expansion of Eloc(l1,l2)"
 811       do i=0,nloctyp-1
 812         read (ifourier,*,end=115,err=115)
 813         read (ifourier,*,end=115,err=115) (b(ii,i),ii=1,13)
 814         if (lprint) then
 815         write (iout,*) 'Type ',onelet(iloctyp(i))
 816         write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
 817         endif
 818         if (i.gt.0) then
 819         b(2,-i)= b(2,i)
 820         b(3,-i)= b(3,i)
 821         b(4,-i)=-b(4,i)
 822         b(5,-i)=-b(5,i)
 823         endif
 824 c        B1(1,i)  = b(3)
 825 c        B1(2,i)  = b(5)
 826 c        B1(1,-i) = b(3)
 827 c        B1(2,-i) = -b(5)
 828 c        b1(1,i)=0.0d0
 829 c        b1(2,i)=0.0d0
 830 c        B1tilde(1,i) = b(3)
 831 c        B1tilde(2,i) =-b(5)
 832 c        B1tilde(1,-i) =-b(3)
 833 c        B1tilde(2,-i) =b(5)
 834 c        b1tilde(1,i)=0.0d0
 835 c        b1tilde(2,i)=0.0d0
 836 c        B2(1,i)  = b(2)
 837 c        B2(2,i)  = b(4)
 838 c        B2(1,-i)  =b(2)
 839 c        B2(2,-i)  =-b(4)
 840 cc        B1tilde(1,i) = b(3,i)
 841 cc        B1tilde(2,i) =-b(5,i)
 842 C        B1tilde(1,-i) =-b(3,i)
 843 C        B1tilde(2,-i) =b(5,i)
 844 cc        b1tilde(1,i)=0.0d0
 845 cc        b1tilde(2,i)=0.0d0
 846 cc        B2(1,i)  = b(2,i)
 847 cc        B2(2,i)  = b(4,i)
 848 C        B2(1,-i)  =b(2,i)
 849 C        B2(2,-i)  =-b(4,i)
 850
 851 c        b2(1,i)=0.0d0
 852 c        b2(2,i)=0.0d0
 853         CCold(1,1,i)= b(7,i)
 854         CCold(2,2,i)=-b(7,i)
 855         CCold(2,1,i)= b(9,i)
 856         CCold(1,2,i)= b(9,i)
 857         CCold(1,1,-i)= b(7,i)
 858         CCold(2,2,-i)=-b(7,i)
 859         CCold(2,1,-i)=-b(9,i)
 860         CCold(1,2,-i)=-b(9,i)
 861 c        CC(1,1,i)=0.0d0
 862 c        CC(2,2,i)=0.0d0
 863 c        CC(2,1,i)=0.0d0
 864 c        CC(1,2,i)=0.0d0
 865 c        Ctilde(1,1,i)= CCold(1,1,i)
 866 c        Ctilde(1,2,i)= CCold(1,2,i)
 867 c        Ctilde(2,1,i)=-CCold(2,1,i)
 868 c        Ctilde(2,2,i)=-CCold(2,2,i)
 869
 870 c        Ctilde(1,1,i)=0.0d0
 871 c        Ctilde(1,2,i)=0.0d0
 872 c        Ctilde(2,1,i)=0.0d0
 873 c        Ctilde(2,2,i)=0.0d0
 874         DDold(1,1,i)= b(6,i)
 875         DDold(2,2,i)=-b(6,i)
 876         DDold(2,1,i)= b(8,i)
 877         DDold(1,2,i)= b(8,i)
 878         DDold(1,1,-i)= b(6,i)
 879         DDold(2,2,-i)=-b(6,i)
 880         DDold(2,1,-i)=-b(8,i)
 881         DDold(1,2,-i)=-b(8,i)
 882 c        DD(1,1,i)=0.0d0
 883 c        DD(2,2,i)=0.0d0
 884 c        DD(2,1,i)=0.0d0
 885 c        DD(1,2,i)=0.0d0
 886 c        Dtilde(1,1,i)= DD(1,1,i)
 887 c        Dtilde(1,2,i)= DD(1,2,i)
 888 c        Dtilde(2,1,i)=-DD(2,1,i)
 889 c        Dtilde(2,2,i)=-DD(2,2,i)
 890
 891 c        Dtilde(1,1,i)=0.0d0
 892 c        Dtilde(1,2,i)=0.0d0
 893 c        Dtilde(2,1,i)=0.0d0
 894 c        Dtilde(2,2,i)=0.0d0
 895         EEold(1,1,i)= b(10,i)+b(11,i)
 896         EEold(2,2,i)=-b(10,i)+b(11,i)
 897         EEold(2,1,i)= b(12,i)-b(13,i)
 898         EEold(1,2,i)= b(12,i)+b(13,i)
 899         EEold(1,1,-i)= b(10,i)+b(11,i)
 900         EEold(2,2,-i)=-b(10,i)+b(11,i)
 901         EEold(2,1,-i)=-b(12,i)+b(13,i)
 902         EEold(1,2,-i)=-b(12,i)-b(13,i)
 903         write(iout,*) "TU DOCHODZE"
 904         print *,"JESTEM"
 905 c        ee(1,1,i)=1.0d0
 906 c        ee(2,2,i)=1.0d0
 907 c        ee(2,1,i)=0.0d0
 908 c        ee(1,2,i)=0.0d0
 909 c        ee(2,1,i)=ee(1,2,i)
 910       enddo
 911       if (lprint) then
 912       write (iout,*)
 913       write (iout,*)
 914      &"Coefficients of the cumulants (independent of valence angles)"
 915       do i=-nloctyp+1,nloctyp-1
 916         write (iout,*) 'Type ',onelet(iloctyp(i))
 917         write (iout,*) 'B1'
 918         write(iout,'(2f10.5)') B(3,i),B(5,i)
 919         write (iout,*) 'B2'
 920         write(iout,'(2f10.5)') B(2,i),B(4,i)
 921         write (iout,*) 'CC'
 922         do j=1,2
 923           write (iout,'(2f10.5)') CCold(j,1,i),CCold(j,2,i)
 924         enddo
 925         write(iout,*) 'DD'
 926         do j=1,2
 927           write (iout,'(2f10.5)') DDold(j,1,i),DDold(j,2,i)
 928         enddo
 929         write(iout,*) 'EE'
 930         do j=1,2
 931           write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i)
 932         enddo
 933       enddo
 934       endif
 935 #endif
 936 C      write (iout,*) 'KURWAKURWA'
 937 #ifdef CRYST_TOR
 938 C
 939 C Read torsional parameters in old format
 940 C
 941       read (itorp,*,end=113,err=113) ntortyp,nterm_old
 942       write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
 943       read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
 944       do i=1,ntortyp
 945         do j=1,ntortyp
 946           read (itorp,'(a)',end=113,err=113)
 947           do k=1,nterm_old
 948             read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
 949           enddo
 950         enddo
 951       enddo
 952       close (itorp)
 953       if (lprint) then
 954         write (iout,'(/a/)') 'Torsional constants:'
 955         do i=1,ntortyp
 956           do j=1,ntortyp
 957             write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
 958             write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
 959           enddo
 960         enddo
 961       endif
 962 #else
 963 C
 964 C Read torsional parameters
 965 C
 966       IF (TOR_MODE.eq.0) THEN
 967
 968       read (itorp,*,end=113,err=113) ntortyp
 969       read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
 970       write (iout,*) 'ntortyp',ntortyp
 971       do iblock=1,2
 972       do i=-ntyp,-1
 973        itortyp(i)=-itortyp(-i)
 974       enddo
 975 c      write (iout,*) 'ntortyp',ntortyp
 976       do i=0,ntortyp-1
 977         do j=-ntortyp+1,ntortyp-1
 978           read (itorp,*,end=113,err=113) nterm(i,j,iblock),
 979      &          nlor(i,j,iblock)
 980           nterm(-i,-j,iblock)=nterm(i,j,iblock)
 981           nlor(-i,-j,iblock)=nlor(i,j,iblock)
 982           v0ij=0.0d0
 983           si=-1.0d0
 984           do k=1,nterm(i,j,iblock)
 985             read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),
 986      &      v2(k,i,j,iblock)
 987             v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
 988             v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
 989             v0ij=v0ij+si*v1(k,i,j,iblock)
 990             si=-si
 991          enddo
 992           do k=1,nlor(i,j,iblock)
 993             read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),
 994      &        vlor2(k,i,j),vlor3(k,i,j)
 995             v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
 996           enddo
 997           v0(i,j,iblock)=v0ij
 998           v0(-i,-j,iblock)=v0ij
 999         enddo
1000       enddo
1001       enddo
1002       close (itorp)
1003       if (lprint) then
1004         write (iout,'(/a/)') 'Torsional constants:'
1005         do i=1,ntortyp
1006           do j=1,ntortyp
1007             write (iout,*) 'ityp',i,' jtyp',j
1008             write (iout,*) 'Fourier constants'
1009             do k=1,nterm(i,j,iblock)
1010               write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
1011      &        v2(k,i,j,iblock)
1012             enddo
1013             write (iout,*) 'Lorenz constants'
1014             do k=1,nlor(i,j,iblock)
1015               write (iout,'(3(1pe15.5))')
1016      &         vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
1017             enddo
1018           enddo
1019         enddo
1020       endif
1021 C
1022 C 6/23/01 Read parameters for double torsionals
1023 C
1024       do iblock=1,2
1025       do i=0,ntortyp-1
1026         do j=-ntortyp+1,ntortyp-1
1027           do k=-ntortyp+1,ntortyp-1
1028             read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
1029 c              write (iout,*) "OK onelett",
1030 c     &         i,j,k,t1,t2,t3
1031
1032             if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
1033      &        .or. t3.ne.toronelet(k)) then
1034               write (iout,*) "Error in double torsional parameter file",
1035      &         i,j,k,t1,t2,t3
1036 #ifdef MPI
1037               call MPI_Finalize(Ierror)
1038 #endif
1039                stop "Error in double torsional parameter file"
1040             endif
1041           read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),
1042      &         ntermd_2(i,j,k,iblock)
1043             ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
1044             ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
1045             read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,
1046      &         ntermd_1(i,j,k,iblock))
1047             read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,
1048      &         ntermd_1(i,j,k,iblock))
1049             read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,
1050      &         ntermd_1(i,j,k,iblock))
1051             read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,
1052      &         ntermd_1(i,j,k,iblock))
1053 C Martix of D parameters for one dimesional foureir series
1054             do l=1,ntermd_1(i,j,k,iblock)
1055              v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
1056              v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
1057              v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
1058              v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
1059 c            write(iout,*) "whcodze" ,
1060 c     & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
1061             enddo
1062             read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),
1063      &         v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
1064      &         v2s(m,l,i,j,k,iblock),
1065      &         m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
1066 C Martix of D parameters for two dimesional fourier series
1067             do l=1,ntermd_2(i,j,k,iblock)
1068              do m=1,l-1
1069              v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
1070              v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
1071              v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
1072              v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
1073              enddo!m
1074             enddo!l
1075           enddo!k
1076         enddo!j
1077       enddo!i
1078       enddo!iblock
1079       if (lprint) then
1080       write (iout,*)
1081       write (iout,*) 'Constants for double torsionals'
1082       do iblock=1,2
1083       do i=0,ntortyp-1
1084         do j=-ntortyp+1,ntortyp-1
1085           do k=-ntortyp+1,ntortyp-1
1086             write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
1087      &        ' nsingle',ntermd_1(i,j,k,iblock),
1088      &        ' ndouble',ntermd_2(i,j,k,iblock)
1089             write (iout,*)
1090             write (iout,*) 'Single angles:'
1091             do l=1,ntermd_1(i,j,k,iblock)
1092               write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
1093      &           v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
1094      &           v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),
1095      &           v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
1096             enddo
1097             write (iout,*)
1098             write (iout,*) 'Pairs of angles:'
1099             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1100             do l=1,ntermd_2(i,j,k,iblock)
1101               write (iout,'(i5,20f10.5)')
1102      &         l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
1103             enddo
1104             write (iout,*)
1105            write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1106             do l=1,ntermd_2(i,j,k,iblock)
1107               write (iout,'(i5,20f10.5)')
1108      &         l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),
1109      &         (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
1110             enddo
1111             write (iout,*)
1112           enddo
1113         enddo
1114       enddo
1115       enddo
1116       endif
1117
1118       ELSE IF (TOR_MODE.eq.1) THEN
1119
1120 C read valence-torsional parameters
1121       read (itorp,*,end=113,err=113) ntortyp
1122       nkcctyp=ntortyp
1123       write (iout,*) "Valence-torsional parameters read in ntortyp",
1124      &   ntortyp
1125       read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1126       write (iout,*) "itortyp_kcc",(itortyp(i),i=1,ntyp)
1127       do i=-ntyp,-1
1128         itortyp(i)=-itortyp(-i)
1129       enddo
1130       do i=-ntortyp+1,ntortyp-1
1131         do j=-ntortyp+1,ntortyp-1
1132 C first we read the cos and sin gamma parameters
1133           read (itorp,'(13x,a)',end=113,err=113) string
1134           write (iout,*) i,j,string
1135           read (itorp,*,end=113,err=113)
1136      &    nterm_kcc(j,i),nterm_kcc_Tb(j,i)
1137 C           read (itorkcc,*,end=121,err=121) nterm_kcc_Tb(j,i)
1138           do k=1,nterm_kcc(j,i)
1139             do l=1,nterm_kcc_Tb(j,i)
1140               do ll=1,nterm_kcc_Tb(j,i)
1141               read (itorp,*,end=113,err=113) ii,jj,kk,
1142      &          v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
1143               enddo
1144             enddo
1145           enddo
1146         enddo
1147       enddo
1148
1149       ELSE
1150 c AL 4/8/16: Calculate coefficients from one-body parameters
1151       ntortyp=nloctyp
1152       do i=-ntyp1,ntyp1
1153         itortyp(i)=itype2loc(i)
1154 c        itortyp(i)=0
1155       enddo
1156       write (iout,*)
1157      &"Val-tor parameters calculated from cumulant coefficients ntortyp"
1158      & ,ntortyp
1159       do i=-ntortyp+1,ntortyp-1
1160         do j=-ntortyp+1,ntortyp-1
1161           nterm_kcc(j,i)=2
1162           nterm_kcc_Tb(j,i)=3
1163           do k=1,nterm_kcc_Tb(j,i)
1164             do l=1,nterm_kcc_Tb(j,i)
1165               v1_kcc(k,l,1,i,j)=bnew1(k,1,i)*bnew2(l,1,j)
1166      &                         +bnew1(k,2,i)*bnew2(l,2,j)
1167               v2_kcc(k,l,1,i,j)=bnew1(k,1,i)*bnew2(l,2,j)
1168      &                         +bnew1(k,2,i)*bnew2(l,1,j)
1169             enddo
1170           enddo
1171           do k=1,nterm_kcc_Tb(j,i)
1172             do l=1,nterm_kcc_Tb(j,i)
1173               v1_kcc(k,l,2,i,j)=-0.25*(ccnew(k,1,i)*ddnew(l,1,j)
1174      &                         -ccnew(k,2,i)*ddnew(l,2,j))
1175               v2_kcc(k,l,2,i,j)=-0.25*(ccnew(k,2,i)*ddnew(l,1,j)
1176      &                         +ccnew(k,1,i)*ddnew(l,2,j))
1177             enddo
1178           enddo
1179 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)
1180         enddo
1181       enddo
1182
1183       ENDIF ! TOR_MODE
1184
1185       if (tor_mode.gt.0 .and. lprint) then
1186 c Print valence-torsional parameters
1187         write (iout,'(a)')
1188      &    "Parameters of the valence-torsional potentials"
1189         do i=-ntortyp+1,ntortyp-1
1190         do j=-ntortyp+1,ntortyp-1
1191         write (iout,'(3a)') "Type ",toronelet(i),toronelet(j)
1192         write (iout,'(3a5,2a15)') "itor","ival","jval","v_kcc","v2_kcc"
1193         do k=1,nterm_kcc(j,i)
1194           do l=1,nterm_kcc_Tb(j,i)
1195             do ll=1,nterm_kcc_Tb(j,i)
1196                write (iout,'(3i5,2f15.4)')
1197      &            k,l-1,ll-1,v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
1198             enddo
1199           enddo
1200         enddo
1201         enddo
1202         enddo
1203       endif
1204
1205 #endif
1206 C Read of Side-chain backbone correlation parameters
1207 C Modified 11 May 2012 by Adasko
1208 CCC
1209 C
1210       read (isccor,*,end=119,err=119) nsccortyp
1211       read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
1212       do i=-ntyp,-1
1213         isccortyp(i)=-isccortyp(-i)
1214       enddo
1215       iscprol=isccortyp(20)
1216 c      write (iout,*) 'ntortyp',ntortyp
1217       maxinter=3
1218 cc maxinter is maximum interaction sites
1219       do l=1,maxinter
1220       do i=1,nsccortyp
1221         do j=1,nsccortyp
1222           read (isccor,*,end=119,err=119)
1223      &nterm_sccor(i,j),nlor_sccor(i,j)
1224 c          write (iout,*) nterm_sccor(i,j)
1225           v0ijsccor=0.0d0
1226           v0ijsccor1=0.0d0
1227           v0ijsccor2=0.0d0
1228           v0ijsccor3=0.0d0
1229           si=-1.0d0
1230           nterm_sccor(-i,j)=nterm_sccor(i,j)
1231           nterm_sccor(-i,-j)=nterm_sccor(i,j)
1232           nterm_sccor(i,-j)=nterm_sccor(i,j)
1233 c          write (iout,*) nterm_sccor(i,j),nterm_sccor(-i,j),
1234 c     &    nterm_sccor(-i,-j),nterm_sccor(i,-j)
1235           do k=1,nterm_sccor(i,j)
1236             read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
1237      &    ,v2sccor(k,l,i,j)
1238             if (j.eq.iscprol) then
1239              if (i.eq.isccortyp(10)) then
1240              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1241              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1242              else
1243              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
1244      &                        +v2sccor(k,l,i,j)*dsqrt(0.75d0)
1245              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
1246      &                        +v1sccor(k,l,i,j)*dsqrt(0.75d0)
1247              v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1248              v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1249              v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1250              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1251              endif
1252             else
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                if (j.eq.isccortyp(10)) then
1258              v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
1259              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
1260                else
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              v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1266              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1267                 endif
1268                endif
1269             endif
1270             v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
1271             v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
1272             v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
1273             v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
1274             si=-si
1275            enddo
1276           do k=1,nlor_sccor(i,j)
1277             read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
1278      &        vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1279             v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
1280      &(1+vlor3sccor(k,i,j)**2)
1281           enddo
1282           v0sccor(l,i,j)=v0ijsccor
1283           v0sccor(l,-i,j)=v0ijsccor1
1284           v0sccor(l,i,-j)=v0ijsccor2
1285           v0sccor(l,-i,-j)=v0ijsccor3
1286           enddo
1287         enddo
1288       enddo
1289       close (isccor)
1290       if (lprint) then
1291         write (iout,'(/a/)') 'Torsional constants of SCCORR:'
1292         do l=1,maxinter
1293         do i=1,nsccortyp
1294           do j=1,nsccortyp
1295             write (iout,*) 'ityp',i,' jtyp',j
1296             write (iout,*) 'Fourier constants'
1297             do k=1,nterm_sccor(i,j)
1298               write (iout,'(2(1pe15.5))')
1299      & v1sccor(k,l,i,j),v2sccor(k,l,i,j)
1300             enddo
1301             write (iout,*) 'Lorenz constants'
1302             do k=1,nlor_sccor(i,j)
1303               write (iout,'(3(1pe15.5))')
1304      &         vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1305             enddo
1306           enddo
1307         enddo
1308         enddo
1309       endif
1310 C
1311 C Read electrostatic-interaction parameters
1312 C
1313       if (lprint) then
1314         write (iout,'(/a)') 'Electrostatic interaction constants:'
1315         write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
1316      &            'IT','JT','APP','BPP','AEL6','AEL3'
1317       endif
1318       read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
1319       read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
1320       read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
1321       read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
1322       close (ielep)
1323       do i=1,2
1324         do j=1,2
1325         rri=rpp(i,j)**6
1326         app (i,j)=epp(i,j)*rri*rri
1327         bpp (i,j)=-2.0D0*epp(i,j)*rri
1328         ael6(i,j)=elpp6(i,j)*4.2D0**6
1329         ael3(i,j)=elpp3(i,j)*4.2D0**3
1330         lprint=.true.
1331         if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
1332      &                    ael6(i,j),ael3(i,j)
1333         enddo
1334       enddo
1335 C
1336 C Read side-chain interaction parameters.
1337 C
1338       read (isidep,*) ipot,expon
1339       if (ipot.lt.1 .or. ipot.gt.6) then
1340         write (iout,'(2a)') 'Error while reading SC interaction',
1341      &               'potential file - unknown potential type.'
1342         stop
1343       endif
1344       expon2=expon/2
1345       write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
1346      & ', exponents are ',expon,2*expon
1347       goto (10,20,30,30,40,50,60) ipot
1348 C----------------------- LJ potential ---------------------------------
1349    10 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
1350       if (lprint) then
1351         write (iout,'(/a/)') 'Parameters of the LJ potential:'
1352         write (iout,'(a/)') 'The epsilon array:'
1353         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1354         write (iout,'(/a)') 'One-body parameters:'
1355         write (iout,'(a,4x,a)') 'residue','sigma'
1356         write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1357       endif
1358       goto 60
1359 C----------------------- LJK potential --------------------------------
1360    20 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
1361      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1362       if (lprint) then
1363         write (iout,'(/a/)') 'Parameters of the LJK potential:'
1364         write (iout,'(a/)') 'The epsilon array:'
1365         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1366         write (iout,'(/a)') 'One-body parameters:'
1367         write (iout,'(a,4x,2a)') 'residue','   sigma  ','    r0    '
1368         write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
1369      &        i=1,ntyp)
1370       endif
1371       goto 60
1372 C---------------------- GB or BP potential -----------------------------
1373    30 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
1374      &  (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip0(i),i=1,ntyp),
1375      &  (alp(i),i=1,ntyp)
1376 C For the GB potential convert sigma'**2 into chi'
1377       if (ipot.eq.4) then
1378         do i=1,ntyp
1379           chip(i)=(chip0(i)-1.0D0)/(chip0(i)+1.0D0)
1380         enddo
1381       endif
1382       if (lprint) then
1383         write (iout,'(/a/)') 'Parameters of the BP potential:'
1384         write (iout,'(a/)') 'The epsilon array:'
1385         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1386         write (iout,'(/a)') 'One-body parameters:'
1387         write (iout,'(a,4x,5a)') 'residue','   sigma  ','s||/s_|_^2',
1388      &       '   chip0  ','    chip  ','    alph  '
1389         write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),sigii(i),
1390      &                     chip0(i),chip(i),alp(i),i=1,ntyp)
1391       endif
1392       goto 60
1393 C--------------------- GBV potential -----------------------------------
1394    40 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
1395      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1396      &  (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1397       if (lprint) then
1398         write (iout,'(/a/)') 'Parameters of the GBV 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  ','    r0    ',
1403      &      's||/s_|_^2','    chip  ','    alph  '
1404         write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1405      &           sigii(i),chip(i),alp(i),i=1,ntyp)
1406       endif
1407       goto 60
1408 C--------------------- Momo potential -----------------------------------
1409    50 continue
1410       read (isidep,*) (icharge(i),i=1,ntyp)
1411 c      write (2,*) "icharge",(icharge(i),i=1,ntyp)
1412       do i=1,ntyp
1413        do j=1,i
1414 c!        write (*,*) "Im in ", i, " ", j
1415         read(isidep,*)
1416      &  eps(i,j),sigma(i,j),chi(i,j),chi(j,i),chipp(i,j),chipp(j,i),
1417      &  (alphasur(k,i,j),k=1,4),sigmap1(i,j),sigmap2(i,j),
1418      &  chis(i,j),chis(j,i),
1419      &  nstate(i,j),(wstate(k,i,j),k=1,4),
1420      &  dhead(1,1,i,j),
1421      &  dhead(1,2,i,j),
1422      &  dhead(2,1,i,j),
1423      &  dhead(2,2,i,j),
1424      &  dtail(1,i,j),dtail(2,i,j),
1425      &  epshead(i,j),sig0head(i,j),
1426      &  rborn(i,j),rborn(j,i),
1427      &  (wqdip(k,i,j),k=1,2),wquad(i,j),
1428      &  alphapol(i,j),alphapol(j,i),
1429      &  (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(i,j),epsintab(i,j)
1430        END DO
1431       END DO
1432 c!      write (*,*) "nstate gly-gly", nstate(10,10)
1433 c! THIS LOOP FILLS PARAMETERS FOR PAIRS OF AA's NOT EXPLICITLY
1434 c! DEFINED IN SCPARM.MOMO. IT DOES SO BY TAKING THEM FROM SYMMETRIC
1435 c! PAIR, FOR EG. IF ARG-HIS IS BLANK, IT WILL TAKE PARAMETERS
1436 c! FROM HIS-ARG.
1437 c!
1438 c! DISABLE IT AT >>YOUR OWN PERIL<<
1439 c!
1440       DO i = 1, ntyp
1441        DO j = i+1, ntyp
1442         eps(i,j) = eps(j,i)
1443         sigma(i,j) = sigma(j,i)
1444         nstate(i,j) = nstate(j,i)
1445         sigmap1(i,j) = sigmap1(j,i)
1446         sigmap2(i,j) = sigmap2(j,i)
1447         sigiso1(i,j) = sigiso1(j,i)
1448         sigiso2(i,j) = sigiso2(j,i)
1449
1450         DO k = 1, 4
1451          alphasur(k,i,j) = alphasur(k,j,i)
1452          wstate(k,i,j) = wstate(k,j,i)
1453          alphiso(k,i,j) = alphiso(k,j,i)
1454         END DO
1455
1456         dhead(2,1,i,j) = dhead(1,1,j,i)
1457         dhead(2,2,i,j) = dhead(1,2,j,i)
1458         dhead(1,1,i,j) = dhead(2,1,j,i)
1459         dhead(1,2,i,j) = dhead(2,2,j,i)
1460         dtail(1,i,j) = dtail(1,j,i)
1461         dtail(2,i,j) = dtail(2,j,i)
1462 c!        DO k = 1, 2
1463 c!         DO l = 1, 2
1464 c!         write (*,*) "dhead(k,l,j,i) = ", dhead(k,l,j,i)
1465 c!         write (*,*) "dhead(k,l,i,j) = ", dhead(k,l,i,j)
1466 c!          dhead(l,k,i,j) = dhead(k,l,j,i)
1467 c!         END DO
1468 c!        END DO
1469
1470         epshead(i,j) = epshead(j,i)
1471         sig0head(i,j) = sig0head(j,i)
1472
1473         DO k = 1, 2
1474          wqdip(k,i,j) = wqdip(k,j,i)
1475         END DO
1476
1477         wquad(i,j) = wquad(j,i)
1478         epsintab(i,j) = epsintab(j,i)
1479
1480        END DO
1481       END DO
1482
1483       if (.not.lprint) goto 70
1484       write (iout,'(a)')
1485      & "Parameters of the new physics-based SC-SC interaction potential"
1486       write (iout,'(/7a)') 'Residues','     epsGB','       rGB',
1487      &  '    chi1GB','    chi2GB','   chip1GB','   chip2GB'
1488       do i=1,ntyp
1489         do j=1,i
1490           write (iout,'(2(a3,1x),1pe10.3,5(0pf10.3))')
1491      &      restyp(i),restyp(j),eps(i,j),sigma(i,j),chi(i,j),chi(j,i),
1492      &       chipp(i,j),chipp(j,i)
1493         enddo
1494       enddo
1495       write (iout,'(/9a)') 'Residues',' alphasur1',' alphasur2',
1496      &  ' alphasur3',' alphasur4','   sigmap1','   sigmap2',
1497      &  '     chis1','     chis2'
1498       do i=1,ntyp
1499         do j=1,i
1500           write (iout,'(2(a3,1x),8(0pf10.3))')
1501      &      restyp(i),restyp(j),(alphasur(k,i,j),k=1,4),
1502      &       sigmap1(i,j),sigmap2(j,i),chis(i,j),chis(j,i)
1503         enddo
1504       enddo
1505       write (iout,'(/14a)') 'Residues',' nst ','  wst1',
1506      &  '    wst2','    wst3','    wst4','   dh11','   dh21',
1507      &  '   dh12','   dh22','    dt1','    dt2','    epsh1',
1508      &  '   sigh'
1509       do i=1,ntyp
1510         do j=1,i
1511           write (iout,'(2(a3,1x),i3,4f8.3,6f7.2,f9.5,f7.2)')
1512      &      restyp(i),restyp(j),nstate(i,j),(wstate(k,i,j),k=1,4),
1513      &       ((dhead(l,k,i,j),l=1,2),k=1,2),dtail(1,i,j),dtail(2,i,j),
1514      &       epshead(i,j),sig0head(i,j)
1515         enddo
1516       enddo
1517       write (iout,'(/12a)') 'Residues',' ch1',' ch2',
1518      &  '    rborn1','    rborn2','    wqdip1','    wqdip2',
1519      &  '     wquad'
1520       do i=1,ntyp
1521         do j=1,i
1522           write (iout,'(2(a3,1x),2i4,5f10.3)')
1523      &      restyp(i),restyp(j),icharge(i),icharge(j),
1524      &      rborn(i,j),rborn(j,i),(wqdip(k,i,j),k=1,2),wquad(i,j)
1525         enddo
1526       enddo
1527       write (iout,'(/12a)') 'Residues',
1528      &  '  alphpol1',
1529      &  '  alphpol2','  alphiso1','   alpiso2',
1530      &  '   alpiso3','   alpiso4','   sigiso1','   sigiso2',
1531      &  '     epsin'
1532       do i=1,ntyp
1533         do j=1,i
1534           write (iout,'(2(a3,1x),11f10.3)')
1535      &      restyp(i),restyp(j),alphapol(i,j),alphapol(j,i),
1536      &      (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(j,i),
1537      &      epsintab(i,j)
1538         enddo
1539       enddo
1540       goto 70
1541 C For the GB potential convert sigma'**2 into chi'
1542       DO i=1,ntyp
1543        chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1544       END DO
1545       IF (lprint) THEN
1546        write (iout,'(/a/)') 'Parameters of the BP potential:'
1547        write (iout,'(a/)') 'The epsilon array:'
1548        CALL printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1549        write (iout,'(/a)') 'One-body parameters:'
1550        write (iout,'(a,4x,4a)') 'residue','   sigma  ','s||/s_|_^2',
1551      &       '    chip  ','    alph  '
1552        write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1553      &                     chip(i),alp(i),i=1,ntyp)
1554       END IF
1555    60 CONTINUE
1556       close (isidep)
1557 C-----------------------------------------------------------------------
1558 C Calculate the "working" parameters of SC interactions.
1559       do i=2,ntyp
1560         do j=1,i-1
1561         eps(i,j)=eps(j,i)
1562           do k=1,4
1563             alphasur(k,j,i)=alphasur(k,i,j)
1564             alphiso(k,j,i)=alphiso(k,i,j)
1565             wstate(k,j,i)=wstate(k,i,j)
1566           enddo
1567           do k=1,2
1568             wqdip(k,j,i)=wqdip(k,i,j)
1569           enddo
1570           do k=1,2
1571             do l=1,2
1572               dhead(l,k,j,i)=dhead(l,k,i,j)
1573             enddo
1574           enddo
1575         enddo
1576       enddo
1577       IF (ipot.LT.6) THEN
1578        do i=1,ntyp
1579         do j=i,ntyp
1580          sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1581          sigma(j,i)=sigma(i,j)
1582          rs0(i,j)=dwa16*sigma(i,j)
1583          rs0(j,i)=rs0(i,j)
1584         enddo
1585        enddo
1586       END IF
1587
1588    70 continue
1589       if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1590      & 'Working parameters of the SC interactions:',
1591      & '     a    ','     b    ','   augm   ','  sigma ','   r0   ',
1592      & '  chi1   ','   chi2   '
1593       do i=1,ntyp
1594         do j=i,ntyp
1595           epsij=eps(i,j)
1596           if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4 .or. ipot.eq.6 )
1597      &    THEN
1598             rrij=sigma(i,j)
1599           else
1600             rrij=rr0(i)+rr0(j)
1601           endif
1602           r0(i,j)=rrij
1603           r0(j,i)=rrij
1604           rrij=rrij**expon
1605           epsij=eps(i,j)
1606           sigeps=dsign(1.0D0,epsij)
1607           epsij=dabs(epsij)
1608           aa(i,j)=epsij*rrij*rrij
1609           bb(i,j)=-sigeps*epsij*rrij
1610           aa(j,i)=aa(i,j)
1611           bb(j,i)=bb(i,j)
1612           IF ((ipot.gt.2).AND.(ipot.LT.6)) THEN
1613             sigt1sq=sigma0(i)**2
1614             sigt2sq=sigma0(j)**2
1615             sigii1=sigii(i)
1616             sigii2=sigii(j)
1617             ratsig1=sigt2sq/sigt1sq
1618             ratsig2=1.0D0/ratsig1
1619             chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1620             if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1621             rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1622           else
1623             rsum_max=sigma(i,j)
1624           endif
1625 c         if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1626             sigmaii(i,j)=rsum_max
1627             sigmaii(j,i)=rsum_max
1628 c         else
1629 c           sigmaii(i,j)=r0(i,j)
1630 c           sigmaii(j,i)=r0(i,j)
1631 c         endif
1632 cd        write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1633           if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1634             r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1635             augm(i,j)=epsij*r_augm**(2*expon)
1636 c           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1637             augm(j,i)=augm(i,j)
1638           else
1639             augm(i,j)=0.0D0
1640             augm(j,i)=0.0D0
1641           endif
1642           if (lprint) then
1643             if (ipot.lt.6) then
1644             write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1645      &      restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1646      &      sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1647             else
1648             write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3),2i3,10f8.4,
1649      &       i3,40f10.4)')
1650      &      restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1651      &      sigma(i,j),r0(i,j),chi(i,j),chi(j,i),
1652      &      icharge(i),icharge(j),chipp(i,j),chipp(j,i),
1653      &     (alphasur(k,i,j),k=1,4),sigmap1(i,j),sigmap2(j,i),
1654      &     chis(i,j),chis(j,i),
1655      &     nstate(i,j),(wstate(k,i,j),k=1,4),
1656      &     ((dhead(l,k,i,j),l=1,2),k=1,2),dtail(1,i,j),dtail(2,i,j),
1657      &     epshead(i,j),sig0head(i,j),
1658      &     rborn(i,j),(wqdip(k,i,j),k=1,2),wquad(i,j),
1659      &     alphapol(i,j),alphapol(j,i),
1660      &     (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(i,j)
1661
1662             endif
1663                 endif
1664         enddo
1665       enddo
1666 C
1667 C Define the SC-p interaction constants
1668 C
1669 #ifdef OLDSCP
1670       do i=1,20
1671 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1672 C helix formation)
1673 c       aad(i,1)=0.3D0*4.0D0**12
1674 C Following line for constants currently implemented
1675 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1676         aad(i,1)=1.5D0*4.0D0**12
1677 c       aad(i,1)=0.17D0*5.6D0**12
1678         aad(i,2)=aad(i,1)
1679 C "Soft" SC-p repulsion
1680         bad(i,1)=0.0D0
1681 C Following line for constants currently implemented
1682 c       aad(i,1)=0.3D0*4.0D0**6
1683 C "Hard" SC-p repulsion
1684         bad(i,1)=3.0D0*4.0D0**6
1685 c       bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1686         bad(i,2)=bad(i,1)
1687 c       aad(i,1)=0.0D0
1688 c       aad(i,2)=0.0D0
1689 c       bad(i,1)=1228.8D0
1690 c       bad(i,2)=1228.8D0
1691       enddo
1692 #else
1693 C
1694 C 8/9/01 Read the SC-p interaction constants from file
1695 C
1696       do i=1,ntyp
1697         read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1698       enddo
1699       do i=1,ntyp
1700         aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1701         aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1702         bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1703         bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1704       enddo
1705
1706       if (lprint) then
1707         write (iout,*) "Parameters of SC-p interactions:"
1708         do i=1,20
1709           write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1710      &     eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1711         enddo
1712       endif
1713 #endif
1714 C
1715 C Define the constants of the disulfide bridge
1716 C
1717 C      ebr=-12.0D0
1718 c
1719 c Old arbitrary potential - commented out.
1720 c
1721 c      dbr= 4.20D0
1722 c      fbr= 3.30D0
1723 c
1724 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1725 c energy surface of diethyl disulfide.
1726 c A. Liwo and U. Kozlowska, 11/24/03
1727 c
1728 C      D0CM = 3.78d0
1729 C      AKCM = 15.1d0
1730 C      AKTH = 11.0d0
1731 C      AKCT = 12.0d0
1732 C      V1SS =-1.08d0
1733 C      V2SS = 7.61d0
1734 C      V3SS = 13.7d0
1735       write (iout,*) dyn_ss,'dyndyn'
1736       if (dyn_ss) then
1737         ss_depth=ebr/wsc-0.25*eps(1,1)
1738 C        write(iout,*) akcm,whpb,wsc,'KURWA'
1739         Ht=Ht/wsc-0.25*eps(1,1)
1740
1741         akcm=akcm*whpb/wsc
1742         akth=akth*whpb/wsc
1743         akct=akct*whpb/wsc
1744         v1ss=v1ss*whpb/wsc
1745         v2ss=v2ss*whpb/wsc
1746         v3ss=v3ss*whpb/wsc
1747       else
1748         ss_depth=ebr/whpb-0.25*eps(1,1)*wsc/whpb
1749       endif
1750
1751 C      if (lprint) then
1752       write (iout,'(/a)') "Disulfide bridge parameters:"
1753       write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1754       write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1755       write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1756       write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1757      & ' v3ss:',v3ss
1758 C      endif
1759       if (shield_mode.gt.0) then
1760       pi=3.141592d0
1761 C VSolvSphere the volume of solving sphere
1762 C      print *,pi,"pi"
1763 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
1764 C there will be no distinction between proline peptide group and normal peptide
1765 C group in case of shielding parameters
1766       VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
1767       VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
1768       write (iout,*) VSolvSphere,VSolvSphere_div
1769 C long axis of side chain
1770       do i=1,ntyp
1771       long_r_sidechain(i)=vbldsc0(1,i)
1772       short_r_sidechain(i)=sigma0(i)
1773       enddo
1774       buff_shield=1.0d0
1775       endif
1776       return
1777   111 write (iout,*) "Error reading bending energy parameters."
1778       goto 999
1779   112 write (iout,*) "Error reading rotamer energy parameters."
1780       goto 999
1781   113 write (iout,*) "Error reading torsional energy parameters."
1782       goto 999
1783   114 write (iout,*) "Error reading double torsional energy parameters."
1784       goto 999
1785   115 write (iout,*)
1786      &  "Error reading cumulant (multibody energy) parameters."
1787       goto 999
1788   116 write (iout,*) "Error reading electrostatic energy parameters."
1789       goto 999
1790  1161 write (iout,*) "Error reading lipid energy parameters."
1791       goto 999
1792   117 write (iout,*) "Error reading side chain interaction parameters."
1793       goto 999
1794   118 write (iout,*) "Error reading SCp interaction parameters."
1795       goto 999
1796   119 write (iout,*) "Error reading SCCOR parameters"
1797       goto 999
1798   121 write (iout,*) "Error reading bond parameters"
1799   999 continue
1800 #ifdef MPI
1801       call MPI_Finalize(Ierror)
1802 #endif
1803       stop
1804       return
1805       end